Characterising hydrothermal fluid pathways beneath Aluto volcano, Main Ethiopian Rift, using shear wave splitting. Journal of Volcanology and Geothermal ,

Geothermal resources are frequently associated with silicic calderas which show evidence of geologically- recent activity. Hence development of geothermal sites requires both an understanding of the hydrothermal system of these volcanoes, as well as the deeper magmatic processes which drive them. Here we use shear wave splitting to investigate the hydrothermal system at the silicic peralkaline volcano Aluto in the Main Ethiopian Rift, which has experienced repeated uplift and subsidence since at least 2004. We make over 370 robust observations of splitting, showing that anisotropy is conﬁned mainly to the top ∼ 3km of the volcanic ediﬁce. We ﬁnd up to 10% shear wave anisotropy (SWA) is present with a maximum centred at the geothermal reservoir. Fast shear wave orientations away from the reservoir align NNE–SSW, parallel to the present-day minimum compressive stress. Orientations on the ediﬁce, however, are rotated NE–SW in a manner we predict from ﬁeld observations of faults at the surface, providing ﬂuid pressures are suﬃcient to hold two fracture sets open. These fracture sets may be due to the repeated deformation experienced at Aluto and initiated in caldera formation. We therefore attribute the observed anisotropy to aligned cracks held open by over-pressurised gas-rich ﬂuids within and above the reservoir. This study demonstrates that shear wave splitting can be used to map the extent and style of fracturing in volcanic hydrothermal systems. It also lends support to the hypothesis that deformation at Aluto arises from variations of ﬂuid pressures in the hydrothermal system. These constraints will be crucial for future characterisation of other volcanic and geothermal systems, in rift systems and elsewhere. © 2018 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).


Background, regional setting and Aluto
Geothermal resources worldwide are very frequently associated with active or recently-active volcanoes (e.g., Glassley, 2010). In these locations, a magmatic heat and fluid source is present, and an accompanying hydrothermal system has developed (e.g., Grant and Bixley, 2011), which permits the extraction of heat for immediate use or electricity generation. For reasons of public benefit and scientific understanding, therefore, it is vital to understand both how such resources may have developed through time as volcanic systems, and the current and future state of the geothermal system. Locations where recently-active volcanoes can be studied in terms of their geothermal potential and volcanic history include the world's major rift systems, such as Iceland (Arnórsson, 1995), northern New White blue-headed arrow shows the present-day extension direction (Saria et al., 2014). Brown lines here and in (c) are faults mapped by Agostini et al. (2011). Panel (c) location indicated by red rectangle. (c) Earthquakes and seismic stations at Aluto. Circles are event locations from Wilks et al. (2017) of earthquakes which yielded at least one splitting measurement of quality '2' or better, or null. Events are coloured by hypocentral depth, where depth is measured below sea level. (Median elevation in this plot is 1924 m above sea level.) Triangles are seismic stations, coloured by the number of good shear wave splitting measurements made at that station. Large black double-headed arrows show the fast orientations of Keir et al. (2011), plotted at the event-station midpoint, with the location of the respective seismic station shown by an inverted black triangle nearby. The black hatched square shows the approximate location of the Aluto-Langano geothermal power plant, which defines the 'centre' of the study region in this work. Black lines are Aluto faults mapped by Hutchison et al. (2015). Note particularly the elliptical caldera fault inferred from fumarole locations. (d) North-south section showing the earthquakes used in this study as circles colour-coded with depth, and the entire catalogue of Wilks et al. (2017) as grey circles. Horizontal axis shows elevation above sea level (asl). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) by fluids moving along recent (<2 Ma) faults which accommodate present-day spreading. These NNE-SSW-striking faults were identified by Mohr (1962) as the Wonji Fault Belt (WFB), and have been identified with the emplacement of dykes. (See Keir et al., 2015 and references therein.) Further signs of Aluto's continuing magmatic activity has come from geodetic observations, which show up to 15 cm vertical displacement over periods of just 6 months centred on Aluto's edifice (Biggs et al., 2011;Hutchison et al., 2016a). In common with other MER volcanoes, this has been occurring since at least 2004, the

ARTICLE IN PRESS
earliest point at which remote observations have been made (Biggs et al., 2011).
Despite the interest in the role of magmatism in rifting in the MER, and the hazards arising from MER volcanoes (e.g., Vye-Brown et al., 2016), little geophysical monitoring is present. This study uses a network of seismometers deployed as part of the Aluto Regional Geophysical ObservationS (ARGOS) project, which included a magnetotelluric (MT) survey (Samrock et al., 2015), global navigation satellite system geodesy, plus geological, petrological and geochemical studies (Hutchison et al., 2015;Gleeson et al., 2017;Hutchison et al., 2016a,b,c;Braddock et al., 2017), alongside the seismic network (Wilks et al., 2017) used here.

Shear wave splitting
Here, we investigate the structure of the volcano-geothermal system beneath Aluto using shear wave splitting in waves from local volcano-tectonic (VT) earthquakes. Splitting occurs when a shear wave travels through an anisotropic medium, and the energy partitions into two orthogonally-polarised waves with different velocities (e.g., Musgrave, 1954). The polarisation of the faster shear wave, 0, and the delay time between the arrival of the two shear waves, dt, can be measured and provide information about the style of anisotropy in the subsurface.
Shear wave splitting is increasingly used to investigate volcanological processes (e.g., Gerst, 2004;Baird et al., 2015, amongst many) often by assuming that temporal changes in the splitting parameters (0 and dt) reflect changes in stress within the volcano. This is on the basis that the subsurface contains microcracks oriented in all directions, and that those microcracks which are oriented favourably open up in the direction of the minimum horizontal compressive stress (e.g., Crampin and Booth, 1985). In this way, 0 will tend to point along the regional maximum horizontal compressive stress, since fractures will be elongated in this direction.
Lithological layering or elongated inclusions (e.g., melt pockets on the metre scale) will also induce splitting if they are smaller than the seismic wavelength. In this case, interpretation of 0 may be more closely related to pre-existing structure and temporal changes may not be visible (Boness and Zoback, 2006).
In either case, the amount of splitting along a given ray path is determined by the density and aspect ratio of cracks or fractures, amongst many other parameters. Because dt depends on the distance travelled by a shear wave, shear wave anisotropy (SWA) is used, which is the amount of splitting normalised by velocity and ray path length.

Previous observations
Anisotropy within the MER has been explored using teleseismic SKS waves (e.g., Ayele et al., 2004;Kendall et al., 2005), but these measurements average out structure over the entire upper mantle and crust. Keir et al. (2011), however, made two measurements of splitting in our study region ( Fig. 1) from local events. The northern observation was made at station E77, located near Lake Ziway, from an event located beneath the lake. The southern observation was made at E79 from an event beneath Lake Langano. They found 0 to be along NNE-SSW directions, using events at 8 km depth, with dt = 0.12 s, corresponding to SWA ≈ 6%. They interpreted their results as due to aligned fractures and faults within dykes, which have been emplaced from the Aluto volcanic centre during the Quaternary period. However, no surface expression of dyking has been recorded in these specific areas (Hutchison et al., 2015;Wolde-Gabriel et al., 1990;Kebede et al., 1984), leaving room for alternative explanations for the observed splitting.
Here we investigate shear wave splitting and fracture-induced anisotropy, for the first time at Aluto, using a much richer dataset from a dense seismic monitoring network.

Data and event locations
The data used in this study are as described by Wilks et al. (2017). They are recordings from 18 three-component Güralp 6TD broadband seismometers located around the edifices of the Aluto and Corbetti volcanoes, with a minimum, mean and maximum station spacing respectively of about 1, 10 and 20 km at Aluto (Fig. 1). The network operated from January 2012 to January 2014. Wilks et al. (2017) located 2162 earthquakes during the network's operation, of which 1361 have been termed 'Aluto events', occurring within 15 km of the centre of the edifice. Events were detected by manual inspection of seismograms and picking of P and S wave arrivals times, which were in turn inverted for the events' onset times and locations. We use these recorded S wave arrivals in this study, which number 1454.
Events were located in two main depth groups: between the surface and sea level (depth of 0 km; 56% of events); and from 4 to 10 km below sea level (bsl; 23%). The remainder (13%) of events are between 0 and 4 km bsl, with a small proportion (7%) deeper than 10 km. Spatially, events cluster around the edifice, which is an effect of both detection bias and true propensity for seismicity to occur along the Artu Jawe fault system. These results are supported by more recent joint hypocentre-velocity inversions (Wilks et al., submitted).

Shear wave splitting analysis
We use the 'minimum-eigenvalue' method of (Silver and Chan, 1991), as modified by Walsh et al. (2013) and implemented in the SHEBA program (Wuestefeld et al., 2010), to retrieve splitting parameters which best linearise the particle motion on the horizontal components of the seismograms. Delay times of the slow shear wave relative to the fast of up to 0.4 s were permitted in the grid search.
To minimise manual processing of the data, we used the multiple window analysis method of Teanby et al. (2004) to automatically choose time windows around the arrival based on 100 trial analysis windows. Initially, start times of windows were trialled between 0.22 and 0.04 s before the S wave onset, with end times of between 0.15 and 0.27 s after the S wave pick. These were determined by manual picking of a subset of the data. All results were manually inspected and classified as quality '1' (the best), '2', '3' or 'null' (clearly no splitting present). In some cases, analysis windows were repicked manually where these automatic times did not properly contain the arriving S wave energy, and the results updated.
Quality '1' results satisfied the following criteria: 1s uncertainty in fast orientation, D0 ≤ 10 • ; 1s uncertainty in delay time D(dt) ≤ 0.01 s; signal-to-noise ratio SNR ≥ 5; clearly elliptical particle motion before correction with optimum splitting operator; clearly linear particle motion thereafter.
Quality '2' results satisfied the same criteria, but with D0 ≤ 20 • , D(dt) ≤ 0.02 s and SNR ≥ 3, whilst for quality '3' results D0 ≤ 30 • and D(dt) ≤ 0.04 s. For the latter, results were only retained where results were unambiguous and clearly positive, but which suffered large uncertainties because of the limited frequency content of the seismic recording in the analysis window (see e.g., Walsh et al., 2013). We recommend that category '3' results be used only in combination with higher-quality results at the same station because of the large stated uncertainties.
Null results were classified when particle motion was clearly linear before analysis, and SNR ≥ 5.

ARTICLE IN PRESS
We limited our analysis to event-station pairs whose straightline angle of incidence was less that 55 • from the vertical. Although this appears to be outside the shear wave window within which shear wave polarisations remain unaffected by interaction with the free surface (e.g., Booth and Crampin, 1985), this angle is an overestimate as we do not perform ray-tracing on each path. P-and S-wave velocities in the near surface are believed to be approximately 3.3 and 1.9 km s −1 respectively (Daly et al., 2008;Wilks et al., 2017), and likely lower still in the top few km, so it is expected that significant deviation of the ray will occur and this range of values will ensure accurate shear wave splitting observations. Manual inspection of a number of seismograms at the largest straight-line incidence angles showed no evidence of a free-surface effect on the shear wave polarisation. This was diagnosed by converting the traces via the free-surface transform (Kennett, 1991) into P, SV and SH components for a variety of assumed subsurface velocities at the predicted ray slownesses and investigating for the removal of any significant P-S coupled energy.

Result quality
Of the 1454 S-waves in the catalogue, 1207 yielded splitting parameters. About 13% were quality '1', 14% quality '2' and the remainder, 73%, '3'. Just 26 results were clearly null. The location of earthquakes giving reliable (qualities '1', '2' and 'null') results and the entire event catalogue (Wilks et al., 2017) are shown in Fig. 1. Fig. 2 shows the distribution of signal-to-noise ratio (SNR) for all acceptable ('1'-'3' and 'null') shear wave splitting measurements made at all stations. We define SNR as where A w is the peak amplitude on the horizontal components within the shear wave splitting analysis window (A = A 2 N + A 2 E 1/2 for the two horizontal components, N and E), and A n is the peak amplitude in the 'noise window' before the S-wave onset. Here we use a noise window of length 1 s. This is a simpler measure than that of Restivo and Helffrich (1999), for whom

Fig. 2.
Observed signal-to-noise ratio (SNR) for all S-waves at all stations (blue bars), and SNR of only those yielding good shear wave splitting measurements (orange), using the SNR measure defined in Section 3.1. SNR is also shown using the measure of Restivo and Helffrich (1999) where R max is the maximum absolute amplitude of the 'R' component within the shear wave splitting analysis window, s T is the standard deviation of the amplitudes on the 'T' component within the window, and both of these quantities are measured on the traces after correction by the best-fitting splitting operator recovered in the analysis. The R component is parallel to the incoming wave polarisation direction at the station, and T is perpendicular thereto. The SNR we use here can be computed without first performing shear wave splitting analysis, which for high-frequency data can be an expensive operation. It also does not depend on the successful retrieval of accurate splitting parameters. We therefore propose our measure as a potential pre-filter before performing analysis. We note that using both SNR and SNR RH , measurements which were manually assigned as 'acceptable' were possible when SNR was at least between 2 and 3 times the mean of the Poisson distribution which fits the observed SNR for all signals. Although outside the scope of this study, this may be of interest for future studies into automatic determination of useful shear wave splitting measurements.

Strength and depth of anisotropy
The average splitting delay time isdt = (0.11 ± 0.06) s. We compute shear wave anisotropy (SWA) according to Thomas and Kendall (2002), assuming straight line raypaths between the event and station, and taking an average velocity along the ray from the 1D model used by Wilks et al. (2017). This leads to values of SWA up to 15%, with a median of 1.2%. (Fig. S1 in the Supplementary information shows the distribution of SWA.) These are likely mild overestimates because the rays will not be straight, though this is less of a problem for the shallower events where less bending will occur. Even for the deepest events with the largest epicentral distances, the overestimation in SWA due to ray bending is only about 1%. Uncertainty in SWA arising from inaccuracy in the velocity model is estimated to be less than 1% based on bootstrap modelling where we randomly perturbed the velocity model at each station separately and recalculated SWA.
There is no significant trend of dt with either event depth or straight-line distance between the event and station ( Fig. 3a-b). SWA also falls off rapidly with event-station distance. Both of these features show that the bulk of the splitting is accrued near the surface, or at least that any coherent anisotropy is present from the surface only to a shallow depth. If coherent splitting occurred throughout the volume sampled by the waves, then SWA would be constant with distance, and dt would increase with distance. The maximum lower limit of anisotropy is likely around 3 km bsl. This is constrained by the fact that, for the events with the steepest incidence angles (≤20 • from the vertical), dt remains approximately constant with event depth, and the shallowest events giving these observations are at ∼3 km bsl. We observe no clear trend in dt or SWA with the ray azimuth or angle relative to the vertical. Under the assumption that all splitting is accrued above 3 km bsl, the mean shear wave anisotropySWA increases to 2.8% for all cases in Table 1. The distribution of SWA values with and without this assumption is shown in Supplementary Fig. S1.
SWA in the upper crust is commonly thought to range between 0 and 5% when caused by cracks (e.g., Crampin, 1987), and the majority (97% of the total) of the values we find lie in this range. However, some of our observations of SWA here (<1%) are over 10%. Although large, similar and indeed larger values have been inferred seismically before, for instance in carbonate hydrocarbon settings (e.g., Potters et al., 1999, who find values up to 20%).
With the supposition that anisotropy is concentrated near the surface, we search for evidence that there are lateral trends in SWA in events shallower than 3 km bsl. By restricting ourselves to the shallower events for this exercise, we are minimising the scatter in SWA arising from paths which only spend a portion of their time within the region responsible for most of the splitting we observe. Fig. 4  shows the variation of SWA for these events, whose median value is 4%. We take a grid of points separated by 250 m in northerly and easterly directions and include all event-station midpoints within a 2 km radius of each grid point. We then average the SWA within that point, rejecting any points where fewer than 3 observation midpoints occur within that bin. We only include raypaths with lengths less than 10 km, and only those defined as 'Aluto' events whose epicentres are within 15 km of the 'centre' of Aluto (defined by the geothermal power plant; Fig. 1). A very similar picture is produced when using grid sizes between 100 m and 1 km, a search radius between 1 km and 5 km, a minimum number of observations per bin between 2 and 5, or raypaths less than 5, 10 or 20 km in length.

ARTICLE IN PRESS
Although caution is needed in interpreting this figure because of uneven azimuthal coverage of ray paths and the tradeoff between depth, path length and SWA, it appears that these events experience stronger splitting beneath the volcanic edifice, up to mean SWA of ∼4%, whilst paths outside the edifice have values of 1 to 2%. The central high-SWA feature is robust to varying the selection criteria across a range of distances, path lengths, incidence angles and grid size as described above. It is not clear if the high-SWA region in the southwest of Fig. 4 is significant, but as events in this region have the longest raypaths on this diagram, it is unlikely to be caused by high apparent SWA due to short path lengths. This region is in fact known for hot springs (Kebede et al., 1984;Hutchison et al., 2015;Braddock et al., 2017), which may suggest an increased contribution to anisotropy from fluids near the surface. A second outlying high-SWA region in the northeast is constrained by only two observations, and we accordingly do not seek to interpret the signal there.  Firstly, the pink circle shows the approximate location and extent of the clay cap atop the geothermal reservoir, inferred from the mineralogy in boreholes (Gianelli and Teklemariam, 1993) and subsurface resistivity from the inversion of magnetotelluric data (Samrock et al., 2015). It is clear that the region of high-SWA encompasses the clay cap and is centred at a similar location.
Secondly, we show in dashed purple lines the −25% contour of V P /V S (=1.32) from the local tomographic model of Wilks et al.
(submitted) at 3 km bsl. There is again a large overlap between the regions defined by this V P /V S contour and that of SWA > 4%.
Finally, we show the E-W-elongated caldera ring fault inferred from volcanic vent locations (Hutchison et al., 2015). In this case, the correlation between high values of SWA and the inferred ring fault is striking.

Fast orientation trend
Fast orientations for the whole dataset are shown in Fig. 5. The overall mean value for fast orientations0 = (13 ± 7) • is not significantly different (at the 95% confidence limit) from either the trend of seismicity (∼15 • ), or the local trend of the Wonji fault belt (∼12 • ) as found by Agostini et al. (2011).0 is also approximately perpendicular to the plate-spreading direction of ∼100 • (e.g., Bendick et al., 2006;DeMets et al., 2010;Saria et al., 2014).
We investigate the lateral variation in 0 (Fig. 6). We show polar histograms for numbered bins of side length 0.4 • where the observations have been grouped laterally by the event-receiver path midpoint. This presentation significantly reduces the scatter by not assuming that the contribution to the splitting comes primarily from the source or receiver location. The bin size is chosen to permit a sufficient number of observations in each bin, such that trends and multimodality can be seen. Bins south of the edifice (numbered 16-30) generally have modal 0 values close to north, whilst those above the edifice and to the north (numbered 1-15) have modal 0 ∼ NE. The difference between these two sets of orientations is significant at the 95% level when tested using the U 2 statistic (Watson, 1961).
The exception to this dichotomy is bin 8. Here, the rays whose midpoints are contained within this bin show a preponderance of fast orientations ∼ENE-WSW. This trend arises primarily from events to the north and northeast of the stations on the edifice, where ray paths traverse the northern part of the edifice. This corresponds to the location of an inferred caldera ring fault mapped by Hutchison et al. (2015), which here strikes approximately ENE-WSW, parallel

Cause of anisotropy
We have shown that shear wave anisotropy of up to 10% in the top ∼10 km beneath Aluto is concentrated in the top few km. This correlates both with the increased rate of seismicity and elevated b-values-up to b = 2.5-above about 1 km below sea level (Wilks et al., 2017). High values of b in the Gutenberg-Richter relationship (Gutenberg and Richter, 1944) imply that there is a preponderance for smaller-magnitude events compared to the CMER regional trend of b = 1.1 (Keir et al., 2006), which in turn suggests that rocks above sea level at Aluto are significantly weaker, or that pore pressures are significantly elevated, or both. Neither of these inferences are surprising given that significant outgassing of CO 2 has been measured at the edifice (Hutchison et al., 2015) related to extensive fumarolic activity (Hutchison et al., 2015;Kebede et al., 1984;Braddock et al., 2017), and the control on this imposed by significant faulting across the volcano (e.g., Kebede et al., 1984;WoldeGabriel et al., 1990;Hutchison et al., 2015). It therefore seems highly likely that, in common with most other crustal settings, the significant seismic anisotropy we see is due to sub-seismic-wavelength cracks in the brittle upper crust (e.g., Crampin and Booth, 1985), held open in this case by supercritical or gas-rich fluids (see Section 4.4). It seems likely here that the regional stress field and the orientation of pre-existing fractures combine to preferentially align the open fractures. The radial scale of the histograms is saturated such that the minimum radius represents a count of five in any bin, and the maximum is not constrained, meaning that the least populous histograms appear intentionally smaller. Dashed line shows inferred caldera ring fault by Hutchison et al. (2015). Supplementary Table S1 gives bin edge coordinates and mean orientations. Notice that the trend shifts from NNE-SSW to NE-SW north of 7.76 • N (bins 1-15), and that significant deviations from the regional trends occur in bin 8. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Off-edifice fast orientations
The fast orientations revealed in this study show that paths which do not sample the edifice itself, to the south and extreme north of the study area, show a NNW-SSE trend (0 = (11 ± 2) • for eventstation midpoints which lie further than 5 km from the centre of the edifice). This orientation is parallel to the lineaments of the 'Wonji Fault Belt' (Mohr, 1962;Fig. 5), a set of smaller-offset faults inside the rift which are thought to have accommodated strain from 2 Ma, especially strain arising from magma emplacement (Ebinger and Casey, 2001). The WFB orientations also reflect the current strain field (Bendick et al., 2006;DeMets et al., 2010), and the dominant trend of T-axes in focal mechanisms of the events used in this study (Wilks et al., 2017). Off-edifice lineaments of fissures and craters mapped at the surface by Hutchison et al. (2015) are also parallel to the WFB.
Although the overall trend as described holds, there are potential variations in the modal off-edifice 0 values. Boxes 24, 18 and 29 may show the influence of the Harorese Rhomboidal Fault System (Corti, 2009, and references therein), where the WFB and older border fault systems are believed to interact. Here, a small number of 0 values are west of north, which coincides with the trend of faults mapped in this region (Acocella, 2007, and Fig. 1).
All these lines of evidence suggest that the fast orientations we observe for paths mostly travelling away from the edifice are caused by the regional stress field allowing cracks to open preferentially

ARTICLE IN PRESS
with their long axes between azimuths of −5 • and 25 • clockwise from north.

On-edifice fast orientations
In contrast to the off-edifice splits, those observations with paths travelling primarily through the volcanic edifice have a NE-SW trend (0 = (29 ± 2) • for paths with midpoints less than 5 km from the centre), parallel to the faults apparent to the east of the study region and called 'border faults' by Agostini et al. (2011) (Fig. 5). These faults are believed to have been primarily active during initial early breakup, at 6-8 Ma in the CMER (e.g., WoldeGabriel et al., 1990) and are oblique to the present-day spreading direction.
The trend of the border faults is not, however, reflected in any major field observations on the edifice itself. Mapping of Aluto reveals that the primary lineaments are still sub-parallel to the WFB, including the Artu Jawe fault zone (AJFZ; Fig. 4) (Kebede et al., 1984;Hutchison et al., 2015Hutchison et al., , 2016c, the major fault within the structure which seems to control the surficial pattern of outgassing, and probably provides the primary pathway for geothermal fluids to circulate from depth. Although there is a lack of evidence for border fault-parallel structures, Hutchison et al. (2015) find that craters and fissures, as well as volcanic vents, have a second modal azimuth of ∼90 • (east-west). This matches the long axis of the inferred caldera ring fault mapped at the surface and inferred from vent locations (Kebede et al., 1984;WoldeGabriel et al., 1990;Hutchison et al., 2015Hutchison et al., , 2016c. Based on this observation, we suggest that the rotation of 0 is due to the interaction of two primary fracture sets in the shallow (<1 km bsl) subsurface.
To show this mechanism might explain the observed rotation of 0, we perform simple modelling of the effect of two vertical fracture sets on the anisotropy of an otherwise isotropic medium, using the theory of Grechka (2007), and the expressions of Hudson (1981) to compute the excess compliances associated with dry (i.e., gas-rich) ellipsoidal cracks. This describes a fractured medium in terms of the dimensionless fracture density of ellipsoidal cracks where N is the number of cracks per volume V, and a is the mean crack radius. For details of the modelling, the reader is referred to Verdon et al. (2009), noting that in this case, because we show values of SWA and 0 for vertical rays, the isotropic background velocities have no effect on our results here. We impose on the background medium one set of fractures parallel to the WFB and to the off-edifice 0, with azimuth 12 • , and a second set with azimuth 90 • . The fracture density of the first, WFBparallel set is fixed at n 1 = 0.15, determined through testing a range of n 1 values, whilst the relative density a of the second set is varied from 0 to 1, giving an absolute fracture density of the second set n 2 = an 1 . a then is a measure of how relatively many E-W cracks are held open in the rock mass through which shear wave splitting is accrued. The results are shown in Fig. 7, which indicates that our model agrees with observations when a is in the approximate range 0.5-0.6 and n 1 = 0.15. (The same is true for the range 0.10 ≤ n 1 ≤ 0.20.) This implies that the secondary E-W fracture set seen at the surface, as well as WFB-parallel fractures, are held open to the depths at which most splitting is accrued, which in this case is < 3 km bsl.
A dimensionless fracture density of 0.15 is high, though by no means unprecedented. For example, recent analysis of fractures in boreholes at a comparable location, the andesitic geothermal field in Rotokawa, New Zealand (Massiot et al., 2017), shows a range of n (P 33 in their notation) up to 0.24. These measurements were made at depths >2 km, similar to the depths at which we believe splitting is accrued at Aluto.

Mechanisms supporting east-west fractures at Aluto
In the locality of Aluto, it is at first glance unintuitive that an east-west fracture set should be created, as the regional stress field, acting alone, clearly produces a minimum horizontal compressive stress subparallel to this. Indeed, we see no evidence of such E-W microstructures away from Aluto's edifice, which generally follow the WFB. Several authors have suggested mechanisms by which edifice loading and magmatic overpressure can create differential stresses which might open up fractures (e.g., Muller and Pollard, 1977;Pinel and Jaupart, 2003;Bagnardi et al., 2013;Muirhead et al., 2015;Wadge et al., 2016), though these generally act either circumferentially or radially.
Alternatively, evidence from sand-box experiments of caldera formation supports the idea that multiple fracture sets necessarily form within the plug of collapsed material (e.g., Walter and Troll, 2001) in order to accommodate the increased space available when ring faults dip outwards, though recent 2D discrete-element method modelling implies inward-dipping faults are also possible (Holohan et al., 2015). Further cycles of deformation, such as observed at Aluto, may enhance these fracture sets further. (See Acocella, 2007 for a review.) Elliptical calderas such as Aluto's have often been interpreted as caused simply by the presence of a differential horizontal stress on collapse (e.g., Bosworth et al., 2003), but it has also been suggested that pre-existing crustal structures determine the long axis (Acocella et al., 2002;Robertson et al., 2016;Wadge et al., 2016). Although there is no clear evidence of this being that case at Aluto, the presence of pre-existing E-W structures would also encourage the creation of the second fracture set that is observed at the surface, if for instance a damage zone is associated with the structure. This is compatible with the observed concentration of anisotropy near the surface, as any E-W fractures would need to be held open by fluids at high pressures at greater depth, which we consider unlikely away from the geothermal system, and so anisotropy from any cross-rift structure would fall away with depth. Pre-existing cross-rift structures are observed to control the hydrothermal system and surface deformation at Corbetti caldera, MER (Lloyd et al., 2018). Fig. 8. Schematic diagram of the hydrothermal and magmatic systems beneath Aluto. The Artu Jawe Fault Zone (AJFZ; grey plane) provides the primary pathway for magmatic fluids (red arrows) to ascend from depth, as well as the route by which most meteoric water (blue arrows) is returned to the surface. Away from the geothermal reservoir (orange spheroid), the AJFZ and the co-planar Wonji faults determine the crack orientations which lead to AJFZ-parallel fast shear wave orientations (shown as large black double-headed arrows beneath red seismic stations). The geothermal reservoir is fed with heat from the magma body (hazy red spheroid) below via the advection of fluids, and contains multiple fracture sets (black lines) held open by gas-rich fluid overpressure, which in turn cause fast shear waves to be oblique to the AJFZ within the caldera. A hypothesised ring fault is shown, delineating the boundary of the fractured geothermal reservoir. The clay cap (green spheroid) insulates the reservoir below, though does not contribute significantly to shear wave splitting. Dashed lines are example seismic ray paths and blue-yellow spheres are earthquakes hypocentres. The surface topography is shown exaggerated three times, but other features are approximately to scale. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

ARTICLE IN PRESS
Even if multiple fracture sets are present, in the presence of the regional stress field alone, only one set would be held open. We suggest therefore that over-pressurised fluids are present beneath the edifice, and the overpressure is sufficient to hold open both fracture sets. This agrees with both the elevated b-values found near the surface (Wilks et al., 2017) and local earthquake tomography showing V p /V S to be in the range 1.45-1.65 (Wilks et al., submitted). Indeed, V P /V S is extremely low (∼1.3) in the region where SWA is largest (Fig. 4). These results further imply that the over-pressurised fluid is gas-rich: V P /V S falls in rocks bearing over-pressurised gas because the bulk modulus is significantly reduced, which reduces V P more than V S (e.g., Ito et al., 1979), and this has been observed to occur in geothermal reservoirs similar to Aluto's (e.g., Lees and Wu, 2000). Hence the combined geophysical observations strongly suggest the presence of gas within the reservoir, which agrees well with the observations of gas-liquid ratios made by Gianelli and Teklemariam (1993). Fig. 8 depicts our suggested mechanism for anisotropy beneath Aluto.

Alternative causes of anisotropy
Whilst upper crustal anisotropy is considered to be most likely caused by aligned, fluid-filled fractures in almost all settings, there are other potential causes of anisotropy which we discuss now.
Firstly, there is ample evidence that the lower crust and upper mantle in the MER contains significant volumes of silicate melt (e.g., Kendall et al., 2005;Hammond et al., 2014), which when aligned in pockets would cause significant anisotropy and high electrical conductivity. However, we rule this out as a cause for the anisotropy we observe, because wells drilled down to 1.5 km bsl at Aluto show no evidence for present-day melt at the relevant depths (Gianelli and Teklemariam, 1993;Gizaw, 1993;Teklemariam, 1996), and on the basis of the MT-derived resistivity structure (Samrock et al., 2015), which shows high resistivity beneath the edifice outside the limited region of an inferred clay cap.
Secondly, the alignment of intrinsically-anisotropic mineral grains would also lead to bulk anisotropy (termed lattice-preferred orientation, LPO). This is a major feature of shale-dominated basins, since shale minerals can be very anisotropic and crystals are typically aligned very well by the depositional process (Valcke et al., 2006;Lonardelli et al., 2007). Mapping of Aluto (Kebede et al., 1984;Hutchison et al., 2015Hutchison et al., , 2016c does not appear to show significant LPO in the dominant eruptive products since much is glassy or pumiceous. LPO may be important, however, where significant alteration of eruptive products has occurred, creating anisotropic clay minerals and potentially aligning the crystals. Teklemariam et al. (1996) inferred a dome of hydrothermally altered clay-rich material beneath Aluto at ∼500 m above sea level (asl) from mineral assemblages in wells drilled across the edifice down to 2.5 km below the surface. Samrock et al. (2015) find a conductive region beneath the centre of the edifice at ∼1 km asl which correlates well with this laterally as well as in depth. Such caps are common in high-enthalpy geothermal systems (e.g., Grant and Bixley, 2011). Primarily, we do not consider the presence of a clay cap itself to be a significant cause of the anisotropy because of its inferred depth-our results do not show a marked increase in SWA at 1 km asl, as would be expected. In addition to this, it is not clear why LPO in this region would produce the observed values of 0. If the fabric had a horizontal foliation, as is likely if the maximum compressive stress is vertical as expected, then all values of 0 would be perpendicular to the backazimuth, which is not the case-the circular mean difference is (23 ± 1) • at 1s. It would also cause a very strong dependence of SWA on incidence angle. If on the other hand the fabric was vertical and had a strike parallel to the AJFZ, then it would only produce fast orientations as for the off-edifice results. In fact, we see a rotation of 0 and we do not observe a strong variation in SWA with ray inclination. Any major contribution to splitting from LPO in the clay cap would have to be explained by a foliation at an angle to the AJFZ, for which we can see no evidence.
Finally, layering of seismically heterogeneous material at a scale smaller than the seismic wavelength would lead to anisotropy. This might arise from the layering of volcanic material. However, lava flows and tephra deposits are laid down subhorizontally, and this

ARTICLE IN PRESS
again would generate shear wave splitting much more in waves travelling near-horizontally, and no splitting in vertical waves, coupled with an SH-fast 0. As before, this is in contrast to what we observe. Frozen igneous intrusions might also be present, and in the case of dykes, would produce values of 0 parallel to their strike. Effective medium modelling (Postma, 1955) requires the edifice to be heavily dyked (over 10% by volume) with dykes of velocities reduced by 10% to match the SWA seen. Dykes with greater velocity than the surrounding medium would generate more splitting, but it is once more not clear why intrusions should strike obliquely to the WFB and current extension direction, rather than radially or circumferentially. Mapping of the surface also does not suggest such extensive intrusion (Kebede et al., 1984;Hutchison et al., 2015Hutchison et al., , 2016c.

Causes of unrest at Aluto
Aluto's pattern of periodic rapid uplift and slower subsidence, since at least 2004 (Biggs et al., 2011;Hutchison et al., 2016a), has been attributed to a number of different causes. Samrock et al. (2015) discuss the possibility that swelling in the clay cap at ∼1 km asl and thermal expansion in the geothermal reservoir at ∼1 km bsl might lead to the observed ground deformation. However, whilst it may be possible that the magnitudes could be reproduced by these mechanisms alone, more recent modelling by Hutchison et al. (2016a) using additional InSAR observations suggests that a source as shallow at 1 km asl could not be responsible, and places the inflationary source at ∼3 km bsl. Similarly, the lateral extent of the resistivity anomaly is not sufficient to reproduce the ground deformation observations. This would argue against any major contribution from the clay cap.
This study suggests that most anisotropy is concentrated above ∼3 km bsl, in a highly fractured, over-pressurised region which likely in part makes up the geothermal reservoir. This correlates with a source of subsidence at 1.5 to 2 km bsl suggested by Hutchison et al. (2016c), whereby the geothermal reservoir deflates over time due to thermal effects, but also the loss of fluids, after the initial inflationary period due to the injection of more volatile-rich silicate melt or magmatic fluid at ∼3 km bsl. This model would predict a temporal variation in the flux of fluids through the near-surface, correlated to the deformation signal, and a concomitant variation in splitting, but regrettably our data are not numerous enough to resolve any temporal trends. We suggest that analysis of a longer time series of shear wave splitting data may be able to address this issue in the future.

Conclusions
Using broadband seismic recordings at Aluto volcano, in the Main Ethiopian Rift, we make approximately 370 high-quality, robust shear wave splitting measurements from local earthquakes up to 40 km deep. We find pervasive splitting which does not increase with depth, showing that shear wave anisotropy of between 0.2% and 10% is present, and is confined to the top 5 km or less. The fast orientations we see outside the main edifice align with the Wonji Fault Belt, a series of faults which accommodate present day strain, consistent with these fractures being used as conduits of fluids in the near surface. Beneath the edifice, observed fast orientations are NE-SW and can be explained by the combination of the two dominant fracture sets observed in the field, one of which is parallel to the current extension direction (E-W), and the other the WFB. Anisotropy varies very strongly laterally, and is highest beneath the edifice, collocated with the geothermal reservoir beneath Aluto which may be bounded by the inferred caldera ring fault. We suggest that overpressure of gas or gas-rich fluids in the geothermal reservoir, sourced from a deeper magma supply, maintains fluid pathways in multiple fracture sets. This study shows the effectiveness of using shear wave splitting as a measure of fracture density, allowing us to image regions of high permeability within the volcano and its hydrothermal system.