The Coupled Magmatic and Hydrothermal Systems of the Restless Aluto Caldera, Ethiopia

Seismicity can be used to better understand interactions between magma bodies, hydrothermal systems and their host rocks—key factors influencing volcanic unrest. Here, we use earthquake data to image, for the first time, the seismic velocity structure beneath Aluto, a deforming volcano in the Main Ethiopian Rift. Traveltime tomography is used to jointly relocate seismicity and image 3D P- and S-wave velocity structures and the ratio between them (V P/V S). At depths of 4–9 km, the seismicity maps the top of a large low velocity zone with high V P/V S, which we interpret as a more ductile and melt-bearing region. A shallow (<3 km) hydrothermal system exhibits low seismic velocities and very low V P/V S (∼1.40), consistent with the presence of gases exsolved from a deeper melt-rich mush body. The Artu Jawe fault and fracture system provides the migration pathway that connects the deeper mush body with the shallow hydrothermal system. Together, these observations demonstrate that the interaction between magmatic and hydrothermal systems, driven by the exchange of fluids, is responsible for the restless behavior of Aluto.


INTRODUCTION
Interactions between magmatic and hydrothermal systems beneath volcanoes are poorly understood, but play a crucial role in eruptions, especially those that are phreatomagmatic (Pritchard et al., 2019;Troise et al., 2019). This interaction is also important for successful geothermal energy production in high-enthalpy volcanic systems (Reinsch et al., 2017). Seismic methods offer insights into the nature of these thermo-magmatic systems (e.g., Greenfield et al., 2016;Hooft et al., 2019;Wespestad et al., 2019). Seismicity indicates fluid movements through faults and conduits (e.g., Prejean et al., 2002;Hudson et al., 2017;Greenfield et al., 2019a), but can also be used to image their subsurface structure. For example, seismic velocities and their ratio (V P /V S ), can be used to image fault structures and regions of partial melt and over-pressured gases (e.g., Johnson and Polland, 2013;Muksin et al., 2013). Such observations have implications for understanding volcanic unrest, assessing volcanic hazard and optimizing geothermal exploration.
The dynamic nature of the volcanoes of the Main Ethiopian Rift (MER) has only recently been realized (Biggs et al., 2011). The population exposure index for these volcanos is high (PEI five or greater), as nearly 10 million people live in the region (Brown et al., 2015). Over the past decade, a number of studies have sought to better understand the dynamics and eruptive histories of these volcanoes. Experiments have revealed a rich complexity of seismicity beneath Corbetti (Lloyd, et al., 2018a;Lavayssière et al., 2019), Bora-Tulu Moye (Greenfield et al., 2019a;Greenfield et al., 2019b) and Aluto (Wilks et al., 2017). Satellite imaging suggests frequently recurring and shallow magmatic activity (Biggs et al., 2011;Hutchison et al., 2016b;Lloyd et al., 2018a), and a role for pre-existing structures in the development of volcanism (Lloyd et al., 2018b). Magnetotellurics (Samrock et al., 2015;Lloyd et al., 2018a;Hübert et al., 2018) and gravity (Gottsmann et al., 2020) offer further insights into melt storage and migration beneath these volcanoes. However, the location of magmatic and hydrothermal reservoirs is still somewhat elusive. Here we focus on Aluto volcano, which lies just south of the town of Ziway, between Lake Ziway and Lake Langano (Figure 1).
Aluto is a silicic volcanic center that is thought to have begun erupting at about 0.5 Ma in association with the development of the regional Wonji Fault Belt (Hutchison et al., 2016a;Agostini et al., 2011), and whose most recent volcanic deposits are dated to 0.4 ka (Hutchison et al., 2016a). Regional earthquake catalogs (Keir et al., 2006) suggest that earthquakes up to magnitude three occur along the relatively low-seismicity Aluto-Gedemsa magmatic segment, but are likely associated with the border faults to the east of our study region ( Figure 1) rather than volcanic activity itself. The volcano has however displayed episodic surface deformation (Biggs et al., 2011;Hutchison et al., 2016b). Two pulses of rapid uplift were observed in 2003 and 2008 (15 cm in 10 months and 10 cm in 6 months), but since then the volcano has been slowly subsiding at rates <3 mm/yr (Biggs et al., 2011;Birhanu et al., 2019).
A series of recent experiments have provided further insights into the active nature of Aluto. Seismicity (Wilks et al., 2017) and observations of fracture-induced seismic anisotropy (Nowacki et al., 2018) reveal a complex fault system and fracture network that extends well below an active hydrothermal system. The pattern of seismicity is seasonal, with the peak of seismicity occurring 2-3 months after the heavy rainy season coincident with lake loading and subsidence (Birhanu et al., 2019). Surface and satellite mapping and CO 2 flux surveys have shown that magmatic and hydrothermal upwelling follows recent faults and fracture networks generated by Quaternary to Recent rifting, namely the Wonji Fault Belt and the associated Artu Jawe Fault Zone (AJFZ) ( ELC Electroconsult, 1986;Hutchison et al., 2015;Hutchison et al., 2016a;Braddock et al., 2017). The high temperature gradient and active hydrothermal system has made Aluto a viable and productive geothermal resource; Ethiopia's first, and at present only, geothermal power plant was built within the caldera in 1999 (Hochstein et al., 2017).
In this study, we image the seismic structure beneath Aluto using local earthquake tomography, providing new insight into the magmatic and hydrothermal processes that drive volcanic unrest. Earthquake locations reveal the dynamic nature of the volcano, whilst seismic velocities can be used to map fluids. Fluidrich rocks exhibit reduced shear moduli, but their bulk moduli are sensitive to the gas content of the fluid. Gases are compressible, whilst liquids are less so. Therefore, tomographic images of V P /V S ratios can be used to further highlight fluid-rich regions, but also provide information about the nature of the fluid. Broad, high V P / V S regions have been identified beneath volcanic arcs (Syracuse et al., 2008), while small pockets of V P /V S >1.8 at depths less than 5 km have been interpreted as individual melt reservoirs beneath volcanoes (Koulakov et al., 2009;Jaxybulatov et al., 2011). Smaller elevated V P /V S anomalies (1.7-1.8), have been interpreted as fluid-rich sediments (Hansen, 2004;Muksin et al., 2013) and shallow steam condensates (Chiarabba and Moretti, 2006). In contrast, shallow, low V P /V S ratios (∼1.6) have been observed within geothermal areas such as Mammoth Mountain and the Geysers, California (Julian et al., 1998;Dawson et al., 2016). Theoretical considerations (Mavko and Mukerji, 1995), experimental results (Dvorkin et al., 1999) and insights from the hydrocarbon industry (e.g., Harris et al., 1996;Yu et al., 2015) show that V P /V S ratios are sensitive to pore-fluid compressibility through its effect on V P, and suggest that low V P /V S regions represent gas-rich reservoirs. Thus, 3-D images of V P /V S are extremely useful in mapping melt (in mush or magma), fluids and gases within geothermal and volcanic environments.

MATERIALS AND METHODS
We use local earthquakes from the Aluto Research and Geophysical ObservationS (ARGOS) event catalog (Wilks et al., 2017) that were recorded at 14 sites between January 2012 and January 2014 (Figure 1), during which time the volcano was slowly subsiding. Event detection was performed manually for P-and S-wave arrivals, with the arrival picks given weightings that reflect the clarity of their onsets (Wilks et al., 2017).
Initially, hypocentral locations and origin times were determined using NONLINLOC (Lomax et al., 2000), a nonlinear probabilistic global search algorithm, as outlined in Wilks et al. (2017). Since seismic tomography is heavily reliant on travel time accuracy, we restrict the number of events we use from this catalog based on NONLINLOC location statistics, pick weightings and a minimum detection threshold of six stations per event; this leaves 161 events and 1,672 arrivals (1,038 P-and 634 S-waves). A 1-D starting velocity model was initially developed by combining well-log data in the uppermost 2-3 km (Gianelli and Teklemariam, 1993;Gizaw, 1993) with a regional tomographic model derived from earthquakes in the Main Ethiopian Rift and Ethiopian Plateau (Daly et al., 2008). This 1-D model is further refined via joint hypocenter and velocity inversion using FMTOMO (Rawlinson and Urvoy, 2006)-see Supplementary Section 2.0, for further details. This refinement reveals that in general the velocity structure beneath Aluto is much slower (up to 20% for P-waves and up to 11% for S-waves; see Supplementary Figure 2.1) than the velocity model for the wider region (i.e., that of Daly et al., 2008).
For the seismic tomography, we use a modified version of the FMTOMO code of Rawlinson and Urvoy (2006), which was originally developed to permit the joint inversion of multiple data types (teleseismic, local earthquake, active source) for several classes of unknowns, including layer velocity, interface geometry and hypocenter location. A grid based eikonal solver known as the Fast Marching Method (Rawlinson and Sambridge, 2003) is used to solve the forward problem of traveltime prediction, and a subspace inversion method (Kennett et al., 1988) is used to adjust the unknowns at each step of the iterative non-linear solution process. Pilia et al. (2013) modified FMTOMO to permit recovery of V P , V S and V P /V S and fully non-linear location of hypocenters (via a grid search). We use this modified version of FMTOMO, but further develop it to produce robust estimates of hypocenter uncertainty based on computing the 95% confidence region around each event. The additional computational effort required to retrieve this information is minimal, because the objective function is already computed at every point during the grid search.

Assessing the Data and Model Fit and Optimization of the Inversion Parameters
For the 3-D seismic tomography we perform six iterations, each of which involves: 1) solution of the forward problem through the current model; 2) inversion for velocity parameters using a 20-D (the number of orthogonal search directions in model space) subspace inversion scheme; 3) hypocenter relocation, based on both P-and S-wave arrival times, using a fully non-linear grid search in the presence of the updated velocity model (Pilia et al., 2013). The inversion for P-and S-wave velocities is done together, along with hypocenter location, since both P-wave and S-wave arrival times are used to constrain velocity and location. The inversion for V P /V S is done separately, since it is a linear inversion performed along the ray paths from the S-wave model.
Structure is defined on a cubic B-spline grid that spans 0.8°h orizontally and 44 km in depth, with an approximately uniform node spacing in all three dimensions of ∼2.2 km (42,527 nodes in total). Within our inversion grid, nodes which are relatively undersampled by the data will take values close to the initial reference model, and appear to have ∼0% deviation. These tend to lie at the edges of the model where few events or receivers are present. This attraction to the reference values is due to two features of our inversion approach: 1) regularization, which penalizes changes to model parameters not required by the data, and 2) the subspace inversion scheme, which only modifies parameters that have an effect on the data fit. For a Frontiers in Earth Science | www.frontiersin.org October 2020 | Volume 8 | Article 579699 more complete description of FMTOMO and its modifications for local earthquake tomography, refer to de Kool et al. (2006), Rawlinson and Urvoy (2006), and Pilia et al. (2013). We consider the trade-off between the complexity of the final model (which is a function of the model roughness and model variance) and the accuracy of the data fitting by performing numerous inversions with a range of damping and smoothing parameters (e.g., Eberhart-Phillips, 1986). This is performed using the following trade-off curves: data fit (variance of the residuals) vs. model roughness, while varying the smoothing parameter; data fit vs. model variance, while varying the damping parameter. The details of how we choose the damping and smoothing parameters are included the Supplementary Material. The final variance reduction for each model is: 56% for the V P model; 49% for the V S model; and 34% for the V P /V S model.

Earthquake Relocation Uncertainty
The non-linear inverse problem that we solve is to find the location of each earthquake such that the differences between the observed and predicted arrival times (through a 3-D velocity model) are minimized. This is formalized by specifying an objective function S(m), where m is a vector that defines the latitude, longitude and depth of the earthquake, such that where d obs is a set of N observed arrival times (from a given hypocenter) with the mean subtracted, C −1 0 is a data covariance or weighting matrix which accounts for data picking uncertainty, and d m (m) is a set of predicted traveltimes with the mean subtracted. By removing the mean, we can eliminate the origin time as an unknown in the inversion, since source location can be constrained by the relative on-set times of a phase at a set of stations.
We locate the minimum of S(m) using a grid-search, which is a fully non-linear technique capable of locating the global minimum of S(m) and providing valuable information on location uncertainty [see Lomax et al. (2000) for a comprehensive review of nonlinear earthquake location methods]. The method works by computing S(m) at a regular grid of points in latitude, longitude and depth. This regular sampling of the objective function in 3-D space means that the point of minimum misfit, which coincides with our preferred earthquake location, can be readily identified. A nested grid search is used to locate the hypocenters. In the first step, a coarse step size of 1 km is used, but once the optimum location is found, the search is repeated in the neighborhood of this point using a step size of 100 m.
One advantage of heavily sampling the objective function across a dense 3-D grid is that hypocenter uncertainty information can be obtained with little additional computational effort. Here, we adopt the approach described by Sambridge and Kennett (1986) for defining the 95% confidence region around the solution m obtained from the grid search. This is achieved by identifying the region of parameter space which satisfies where N is the number of data picks and χ 2 3 (0.95) 7.815 is the Chi-square statistic for a p-value of 0.05 and three degrees of freedom. This approach is valid since the objective function 1, with [C −1 0 ] ij δ ij σ ij (where δ ij 1 when i j and δ ij 0 otherwise), can be written as which has a chi-square distribution with N-M degrees of freedom, where M 3 is the number of model parameters. Equation 3 accounts for the effect of errors in the velocity model (which are usually unclear) by rescaling S so that S( m) N − 3, which is simply the expectation value of a chisquare distribution with N−3 degrees of freedom. This is obtained by multiplication with the factor (N − 3)/S( m). Both the objective function 1 and method for computing confidence intervals assume that the data errors are Gaussian, which will influence the final hypocenter locations and uncertainty estimates.
In our application of this approach, we use the boundary of the 95% confidence region as a proxy for location uncertainty and summarize our results in terms of uncertainty in depth, latitude and longitude, as shown in Figure 2. As expected, uncertainty in depth tends to be larger in comparison to uncertainty in latitude and longitude, although for most events it is well below 5 km. A number of events cluster near the upper boundary of the model, and of these a noticeable proportion have depth uncertainties that are larger than average. This occurs because shallow events that are not well constrained in depth cannot be relocated above the model, and hence are forced to reside immediately below the surface. The upper surface of the model is set to be 2.5 km above sea level, which means that all receivers and sources lie within the model volume. While it would be possible to exclude these events from the inversion, they are no less valid than deeper events, which have a similar location uncertainty. Overall, these results suggest that the distribution of seismicity that we interpret in terms of volcanic processes is robust, including the observation that the zone of earthquake activity tends to deepen as we move away from the volcanic center.

Checkerboard Resolution Tests
Since geophysical inverse problems are almost invariably underdetermined and hence solutions are nonunique, it is useful to evaluate the robustness of model solutions derived by seismic tomography. Inherent uncertainties in the data and varying raypath coverage mean that it is important to assess which features of a model are required by the data and which are not, so that a model's robustness can be quantified. Therefore, using the optimized regularization parameters described previously, we evaluate solution robustness using the so-called "checkerboard test" (Glahn and Granet, 1993;Rawlinson and Sambridge, 2003). In this test, the source-receiver configurations of the observational dataset are used to compute a synthetic travel Frontiers in Earth Science | www.frontiersin.org October 2020 | Volume 8 | Article 579699 time dataset. This is calculated through an artificial "checkerboard" structure, which is overlain on the original input model. Gaussian noise with a standard deviation of 20 ms is also added to the synthetic datasets to simulate the assumed picking errors in the "real" data. The checkerboard comprises a 3-D region of alternating high and low velocity anomalies, where the perturbations are set to be ±0.5 km/s for V P (±6.6-16.2%) and ±0.2 km/s (±4.8-9.6%) for V S . The V P /V S input checkerboard is obtained by dividing the V P checkerboard by the V S checkerboard. The tomography is then run with the 1-D reference model as the starting model (whilst also relocating hypocenters) in an attempt to recover the checkerboard pattern. Regions of the solution model that sufficiently recover the checkerboard pattern are considered to be "well-resolved." A known weakness of checkerboard tests is that it is possible for small-scale structures to be well retrieved in comparison to larger scale structures (Léveque et al., 1993); consequently, we generate checkerboards with differing scale lengths of ∼8.8, ∼4.4, and ∼2.2 km. By adopting this approach, we are well positioned to interpret both larger and smaller scale structure in regions of good recovery. We also include a spike test to determine the resolution of discrete velocity anomalies. The input spike is comprised of a low V P , high V S and low V P /V S anomaly centered at approximately 7.78N, 28.8E and ∼1 km depth. The input checkerboard model for V P /V S is shown in Figure 3 and the recovered model is shown in Figure 4. Details of the synthetic recovery tests for V P and V S velocity structure, and the spike test, are presented in Supplementary Section 4 ( Supplementary Figures 4.1-4.5). The results of the synthetic testing using a checkerboard structure for the 161 events is presented in horizontal and vertical slices. For the coarse checkerboard ( Figure 4A), it is evident that V P /V S anomalies are well resolved within the seismic network and down to at least 10 km depth. Outside of the array and below 10 km depth, however, smearing distorts the recovery and highlights that these regions of the model are much less well resolved. The same features are evident in the intermediate checkerboards ( Figure 4B), where the lack of events and poor raypath coverage restrict resolution outside of the network and deeper than 8 km. At a finer scale [ Figure 4C], the checkerboard is even less well defined with restricted recovery below sea level and toward the periphery of the array. Cumulatively this suggests we cannot recover structure on length scales finer than 3 km. With respect to the spike test, as expected, the recovered spike exhibits a decrease in amplitude and a larger footprint compared to the input spike (Supplementary Material 4.5); in the E-W direction this is minimal, but is more noticeable in the N-S direction, which is consistent with the raypath coverage. Nevertheless, this test clearly shows that the anomalies recovered at this location in the final model are likely correct.
The observed lack of resolution at depth is caused by the predominance of events occurring between the surface and 10 km depth, with the region beneath relatively aseismic in comparison. Lateral smearing is particularly evident in the north-south depth slices due to the high concentrations of events to the north and south of the caldera and the increased density of ray paths orientated in that direction. The poor recovery of the Frontiers in Earth Science | www.frontiersin.org October 2020 | Volume 8 | Article 579699 checkerboard at the edges of the east-west cross-sections is caused by the relative lack of events to the east and west of the network (i.e., toward the border faults). At all scale lengths, the amplitudes of the perturbations are underestimated in comparison to the input checkerboards. This occurs due to the damping and smoothing regularization parameters implemented in the inversion, which favor a conservative solution (in terms of amplitude and wavelength of structure) unless required by the data.

RESULTS
Over a 2-year period between 2012 and 2014, Wilks et al. (2017) located 1,361 earthquakes within 15 km of the center of the caldera. As discussed, in an effort to use only the best located events for imaging and to stabilize the inversions, we use a subset of these data comprised of 161 events and 1,672 source-receiver raypaths. Nevertheless, the resulting patterns in seismicity are similar between the original 1,361 events and the final relocated 161 events. Figure 5 shows how the relocated events change with respect to the original locations and Figure 2 shows the uncertainties in the new event locations. In agreement with Wilks et al. (2017), most events are located above sea level, beneath the volcano edifice. Beneath this lies a gap in seismicity until a depth of roughly 5 km. Figures 2 and 5 show how this seismicity deepens to the north and south, capping what Wilks et al. (2017) interpret as a ductile deeper region. There is a weak bimodal distribution, with more events lying beneath the north and south rim of the volcano. Furthermore, the events cluster along the north-south trend of the AJFZ, which was a particularly clear feature in the original seismicity catalog. Figure 6 shows the V P /V S model obtained from the inversion; the V P and V S models are presented in the Supplementary Material, but Figure 2 shows two cross-sections through the V P model. Poorly constrained regions of the model are masked according to the sum of the absolute values of the Fréchet derivatives, which is a similar approach to the so-called Derivative Weight Sum method (e.g., Biryol et al., 2013). Reduced seismic velocities are more prominent in the P-wave model and in the shallow regions. Low V P /V S ratios of 1.45-1.65 (a Poisson's ratio of ] 0.05-0.21) and a negative V P anomaly (up to −5%) extend across the entire range of the volcano that lies above sea level (Figures 2 and 6). The checkerboard tests (Figures 3 and 4; Supplementary Section 3) show that this is a Frontiers in Earth Science | www.frontiersin.org October 2020 | Volume 8 | Article 579699 robust and well resolved feature. Deeper than this, average to elevated V P /V S ratio are observed. It is particularly high (>2.0) in the well resolved region beneath the central stations of the caldera between 5 and 8 km deep ( Figure 6). Furthermore, we observe an additional high V P /V S region to the south of the caldera at ∼3 km depth ( Figure 6). We note that this region is located toward the periphery of the receiver array, but our synthetic resolution tests suggest that it is reasonably well resolved. Due to a lack of deep seismicity, resolution in these models decreases rapidly below ∼10 km. In summary, a shallow region above sea level and beneath the volcano is characterized by abundant seismicity that follows the trend of faults and fractures in the region. It exhibits low seismic velocities and a low V P /V S ratio. Beneath this lies a region of decreased seismicity with evidence for a narrow high-V P /V S region beneath the volcano. There is also evidence for a high V P / V S region to the south of the volcano. Deeper seismicity marks a brittle region that envelops a more ductile region, deepening to south and north of the caldera. There is some evidence for high V P /V S in this region, but resolution is lacking. In the next section we interpret the shallow region in terms of an active hydrothermal system and the deeper region as magmatic in nature.

DISCUSSION AND INTERPRETATION
A conceptual interpretation of the results is shown in Figure 7.
Here we discuss characteristics of the hydrothermal system and the deeper magmatic system and how they interact.

Hydrothermal System
We interpret the region of low V P /V S ratios that lies above sea level in terms of an over-pressurized gas-rich volume (e.g., Dvorkin et al., 1999;Lees and Wu, 2000) occupying highly fractured and hydrothermally altered volcanic products. Low Poisson's ratios are likely amplified by the presence of aligned cracks, fractures and faults, as Poisson's ratios below 0.1 are uncommon in unfractured isotropic materials (Walsh, 1965). Such structural heterogeneities have been shown to decrease Poisson's ratio and V P /V S in low stress environments. We suggest that structures such as ring faults and the AJFZ contribute to the low ratios that we observe (Hutchison et al., 2015). These are linked to fumarolic activity around the volcanic edifice and at the surface (see Figures 1; Braddock et al., 2017). This inference is also supported by shearwave splitting observations, which imply that the volcanic center above sea level is highly fractured due to at least two dominant fracture sets (Nowacki et al., 2018). The sharp transition from low FIGURE 4 | Checkerboard recovery results for V P /V S of differing checkerboard scale lengths (A-C). Compare with input checkerboard structure shown in Figure 3. See Supplementary Figures 4.1-4.5 for equivalent V P and V S checkerboard test results.
Frontiers in Earth Science | www.frontiersin.org October 2020 | Volume 8 | Article 579699 to higher V P /V S at around sea level may also be influenced by the change in lithology from the silicic volcanic rhyolites of Aluto's recent effusive volcanism and lacustrine sediments to the underlying unit of Bofa Basalt, as observed in well log data (Teklemariam, 1996;Teklemariam et al., 1996). The magnetotelluric (MT) study of Samrock et al. (2015) shows a high conductivity (1-2 Ω m) layer above sea level, which they interpret as a hydrothermally altered clay cap. Experiments show that the compressibility of supercritical fluids bears a greater resemblance to steam than liquid and can also produce low V P /V S anomalies in the subsurface (Burnham et al., 1969). At Aluto, deep well data (to 2,500 m below the surface) have shown that temperatures reach a maximum of 360°C at pressures greater than 22 MPa (Gizaw, 1993), conditions at or above the threshold required for water to become supercritical (e.g., Chouet et al., 2008;Zollo et al., 2008). Under the assumption that this maximum temperature is exceeded, we suggest that supercritical fluids might also contribute to the low V P /V S anomaly.

Magmatic System
High V P /V S ratios in volcanic regions are commonly associated with the presence of partial melt in the crust (e.g., Hammond and Kendall, 2016). Fluids have a shear modulus of zero, which reduces V S in partially molten rocks. In the absence of gas there is little change in fluid compressibility and therefore FIGURE 5 | Map of original and relocated hypocenters found by joint inversion for velocity structure and earthquake location. Small white circles show the starting location from Wilks et al. (2017), whilst larger circles show the final relocations. Lines connect the original and relocated hypocenters. Color shows the depth of the relocations. The yellow star shows the inferred subsidence source found by Hutchison et al. (2016b). Red dashed lines indicate the location of the north-south and eastwest profiles shown in Figure 2.
Frontiers in Earth Science | www.frontiersin.org October 2020 | Volume 8 | Article 579699 little change in V P . Hence, the high V P /V S regions beneath sea level in our inversion results suggests a region of partial melt at depth beneath the hydrothermal system. A crystalline mush zone has also been suggested by Gleesen et al. (2017), who through phase equilibria modeling proposed a magma storage depth of 5.6 ± 1 km below the surface. However, it is difficult to convert V P /V S ratios to melt volume fraction, as seismic velocities are as sensitive to the shape of melt inclusions as they are to volume fraction (Hammond and Kendall, 2016). Due to the underlying physics and differences in array geometries, there are often differences between seismic and magnetotelluric images of melt reservoirs (e.g., Whaler and Hautot, 2016;Pritchard et al., 2018). Seismic velocities are sensitive to volume fraction of melt and the shape of the melt inclusions (wetting angle) whereas MT is most sensitive to the connectivity, composition and amount of melt. Previous MT studies show the Aluto upper crust that lies beneath sea level to be relatively resistive (>100 Ω m) and hence melt-poor (Samrock et al., 2015;Hübert, et al., 2018), whereas we see localized regions of high (>2.0) V P /V S values suggesting the presence of partial melt. A suggested explanation for some discrepancy between MT and seismic results is that the melt exists in isolated pockets which are not sufficiently connected for the rock to be conductive. Alternatively, low temperature peralkaline melts with a low wt% H 2 O are relatively resistive and consistent with the observed values (Guo et al., 2016;Hübert, et al., 2018). Such low water content could be achieved following repeated cycles of melting and recrystallisation within a long-lived crystal mush zone. However, the more recent study of Samrock et al. (2020), which reprocessed the same data using the method of Grayver et al. (2019), but accounts for topography and galvanic distortion, shows evidence of moderately conductive (10-20 Ω m) regions in an otherwise resistive crust below sea level. These values can be explained with higher temperature and melt fractions.
The hypocentral relocations improve the delineation of the convex lower boundary of seismicity, which corresponds to the brittle-ductile transition. This boundary is roughly 5 km beneath the edifice, but deepens to 15 km to the north and south ( Figure 7G). Our error analysis of the relocations suggest, with most events having <3 km uncertainties both laterally and in depth, that the pattern of seismicity is robust. While a contribution to hypocenter depth uncertainty may also come from the velocitydepth trade-off, the simple convex shape of seismicity spans the entire model in the N-S direction, which would require a very broad, coherent and high amplitude velocity anomaly to explain the convexity (noting that at 10 km depth, a 1 km increase in hypocenter depth would require a ∼10% increase in the overlying velocity, assuming a background velocity of 5 km/s). As such, we propose that the general increase in depth of the seismicity zone away from the volcanic center is meaningful. Deepening of seismicity away from volcanic centers has also been observed at Corbetti (Lavayssière, et al., 2019) and at Fentale and Dofen, which are volcanoes in the northern MER (Keir et al., 2006), where earthquakes occurred between 8 and 10 km depth beneath each volcano but as deep as 16 km elsewhere in the magmatic segment. Similarly, Hudson et al. (2017) see a similar pattern directly beneath Bárðarbunga caldera in Iceland. A lack of seismicity more than 10 km below the center of the volcano implies a hot, ductile crust, which in turn may suggest a melt-rich mush storage region, consistent with petrological observations (Hutchison, et al., 2016a). There is also some suggestion of a conductive body at these depths beneath Aluto (Samrock et al., 2020). We also see some evidence for high V P /V S ratios in this region, but note that the resolution of our V P /V S inversions is poor below 10 km depth. However, these results together with the convex pattern of deep seismicity observed at Aluto suggests that the seismicity delineates the upper boundary of a region of higher heat flux (spanning ∼20 km laterally), which has a more ductile rheology and little seismicity ( Figure 7H). This is also consistent with the models of the upper crustal plumbing systems at other rift calderas such as Corbetti (Gottsmann et al., 2020).

Coupled Magmatic and Hydrothermal Systems and Uplift
Aluto is a rapidly deforming volcano which shows episodic ground deformation (Biggs, et al., 2011;Hutchison et al., 2016b). Hutchison et al. (2016b) proposed that the uplift was caused by injection of new material into the magmatic system, and subsidence was associated with the magmatic degassing and depressurization of the hydrothermal system. Based on MT results that show a clay cap ∼500 m thick and ∼2-3 km wide, Samrock et al. (2015) also propose that swelling clays and freshwater incursion may be agents for uplift and subsidence at Aluto. Our model, however, suggests the presence of velocity anomalies of significantly larger dimensions and hence supports the hypothesis that the driver of deformation is deeper in origin. Proximal seismicity and high V P /V S ratios delineate a narrow region that connects a deeper magmatic body to the hydrothermal system. Support for this comes from carbon isotope sampling, which shows that Aluto's magmatic and hydrothermal systems are physically connected, where deep (>2 km), hot (>250°C) geothermal fluids receive ongoing input of magmatic volatiles from beneath (Hutchison et al., 2016b). We therefore propose that it is episodic pulses of magma injection at depth that drives this volatile release, which then causes the uplift-subsidence events observed at the surface in this coupled system. Similar mechanisms have been proposed at other volcanoes (e.g., Bárðarbunga, Iceland- Hudson et al., 2017;Corbetti, Ethiopia-Lloyd et al., 2018b).
Reservoir seismicity occurs with deformation, which can be associated with inflation (e.g., Stork et al., 2015) or deflation (e.g., Segall, 1989;van Thienen-Visser and Breunese, 2015), both of which can be at play in a coupled magmatic and hydrothermal system. A sudden release of volatiles from the deeper magmatic system can lead to over-pressured gases within the hydrothermal system (Battaglia et al., 2006), which in turn leads to inflation. The subsequent outflow of volatiles and hydrothermal fluids from the geothermal system generates subsidence at the surface (De Natale FIGURE 7 | Schematic diagram of the interpreted subsurface structure and processes at Aluto, based on the new results and constraints from previously published studies. The key features are: (A) high seismicity within a low V P /V S region of over-pressurized gas at/close to a supercritical state, associated with (B) fumarolic activity and causing (C) outflow driven surface deflation; (D) a high V P /V S region with a significant melt component that has ascended from below; (E) shallow high V P /V S regions of steam condensates that release volatiles to the surface via (F) hot springs; (G) a deeper seismogenic region, convex with latitude; (H) a ductile reservoir of magmatic mush.
Frontiers in Earth Science | www.frontiersin.org October 2020 | Volume 8 | Article 579699 et al., 2001), which at Aluto was observed for the majority of the seismic experiment (Wilks et al., 2017). The sources responsible for significant (>4.3 cm) subsidence prior to 2012 have been located in the uppermost 2-3 km (Biggs et al., 2011;Hutchison et al., 2016b), which correlates with the base of our low V P /V S region ( Figure 6). We note that there may also be subsidence effects associated with cooling. Campi Flegrei, another young, geothermally-active caldera in a state of unrest, exhibits similar gaseous controls in a shallow seismogenic zone (Di Vito et al., 1999;Chiarabba and Moretti, 2006;Battaglia et al., 2008;Chiodini et al., 2012;De Siena et al., 2017). The values of the low V P /V S anomaly are similar in magnitude to those observed at Aluto (∼1.45), with an accompanying low V P anomaly, estimated to be the primary source location of subsidence. At Campi Flegrei, it has been suggested that supercritical fluids beneath a more rigid clay cap are responsible for uplift (Vanorio and Kanitpanyacharoen, 2015). We suggest that shallow hydrothermal processes at Aluto are comparable to Campi Flegrei in terms of the structure and the processes that drive seismicity and deformation. However, the low V P /V S anomaly at Aluto seems to be within or above the clay cap (Samrock et al., 2015). We attribute high-pressure gas phases and the subsequent outflow from the geothermal reservoir as the primary cause of subsidence in the absence of any eruptive behavior. As these are the first seismic images of the volcanoes of the MER, it is not clear how unique Aluto is with respect to other volcanoes in the MER.
The high V P /V S ( Figure 6) region to the south of the caldera at 3 km depth may represent steam condensates ( Figure 7E) that manifest at shallower depths, away from the main volcanic edifice where temperatures are reduced (Aster and Meyer, 1988;Simiyu, 1999). Samrock et al. (2020) observe a highly conductive anomaly at the same depth but more to the east of the volcano. Condensates may form brines that then migrate toward the surface due to increasing pore pressure along fracture networks associated with the AJFZ (Hutchison et al., 2016b). These may feed or heat the hot springs seen at the surface near Lake Langano (Figures 1 and 7F) (Kebede et al., 1985;Hutchison et al., 2015;Braddock et al., 2017;Hochstein et al., 2017). This is also compatible with the elevation-driven, southerly groundwater flow direction from Lake Ziway to Lake Langano (Bernacsek et al., 1992).
Aluto, one of numerous calderas in the MER, is a surface expression of magma-assisted continental rifting. As rifting progresses, strain is localised in magmatic segments rather than the border faults. However, how magma at depth (e.g., as tomographically imaged, for example, by Bastow et al., 2005 andGallacher et al., 2016) feeds the axial magmatic segments of the MER is still to be determined.

CONCLUSIONS
Our study has produced the first seismic images of Aluto, a geothermally active volcano in the Main Ethiopian Rift, which was only recently discovered to be in a state of unrest. These images, coupled with the relocated earthquake hypocenters, provide constraints on the plumbing system beneath the volcano, which has provided new insight into its possible preeruptive behavior. Major findings include: 1) a volatile-rich ductile "magmatic mush" region below 10 km depth, which is capped by a convex layer of seismicity; 2) a region of elevated V P / V S at intermediate depths (5-8 km), also defined by an absence of seismicity, which likely represents the storage of partially molten material from below; 3) the presence of shallow and localized high V P /V S zones away from the volcano, which may represent concentrations of steam condensates that release volatiles to the surface; and 4) a shallow volume of rock (above sea level) containing over-pressurized gas at or close to supercritical conditions (defined by abundant seismicity, low V P and low V P /V S ). Together, these findings demonstrate that unrest at Aluto is driven by the coupling between magma ascent and hydrothermal response, which produce the inflationary and deflationary episodes that have been recently observed. This model is likely applicable to other hydrothermally active volcanoes in a state of unrest.

DATA AVAILABILITY STATEMENT
The seismic dataset analyzed for this study is open access and available on IRIS (Network Code XM). See Wilks et al. (2017) for more information on the ARGOS seismic project. Files containing the final V P , V S and V P /V S solution models and the relocated hypocenters are available at https://doi.org/10.6084/m9. figshare.5892775.v4.

AUTHOR CONTRIBUTIONS
MW processed the data, performed the analysis and wrote the first draft of the paper, as part of his PhD, which was supervised by JK and JW. JK and JB formulated the project. NR developed the inversion code and provided guidance on its use and implementation, and interpretations. MW, JK, AA, AN, and JB acquired the field data. All authors contributed to the interpretation of the results and the writing of the paper. MW, NR, and AN prepared the Figures.

FUNDING
The Bristol University Microseismic ProjectS (BUMPS) provided funding for the seismic experiment and the equipment was loaned from SEIS-UK, a NERC Service and Facility (GEF loan 962