Ductility and Compressibility Accommodate High Magma Flux Beneath a Silicic Continental Rift Caldera: Insights From Corbetti Caldera (Ethiopia)

Large silicic magma reservoirs preferentially form in the upper crust of extensional continental environments. However, our quantitative understanding of the link between mantle magmatism, silicic reservoirs, and surface deformation during rifting is very limited. Here, we focus on Corbetti, a peralkaline caldera in the densely populated Main Ethiopian Rift, which lies above a focused zone of upper mantle partial melt and has been steadily uplifting at a maximum rate of 6.6±1.2 cm yr−1 for more than 10 yr. Numerical modeling shows that a maximum concomitant residual gravity increase of 9±3 μGal yr−1 by the intrusion of mafic magma at ∼7 km depth into a compressible and inelastic crystal mush best explains the observed deformation and gravity changes. The derived magma mass flux of ∼1011 kg yr−1 is anomalously high and at least 1 order of magnitude greater than the mean long‐term mass eruption rate. This study demonstrates that periodic and high‐rate magmatic rejuvenation of upper‐crustal mush is a significant and rapid contributor to mature continental rifting.


Introduction
Magma generation and transport is widely acknowledged as a key process during continental rifting and lithospheric breakup (Ebinger et al., 2017;Karakas & Dufek, 2015). In particular, analysis of the forces required for extension shows the importance of the thermomechanical effects of intrusions in weakening the lithosphere (Buck, 2006). However, the current two-dimensional models consider basaltic dyke intrusions, using parameters based on observations from, for example, Afar and Iceland, where the crust is already thermomechanically weakened and (at least partially) oceanic (Bialas et al., 2010). The patterns of magma storage and transport in rifts where the crust is still dominantly continental are noticeably different (Abebe et al., 2007;Ebinger et al., 2017) but not yet sufficiently well quantified for inclusion in numerical models. In particular, the role of magma and other fluids during continental breakup and the variability in the distribution of magmatism and volcanism along the rift axis (Ebinger et al., 2017) need to be considered.
The Eastern Branch of the East African Rift System (EARS) is an ideal natural laboratory for multiparameteric studies of the magma pathways associated with silicic reservoirs, because (1) 20-30% of the crust comprises new magmatic material (Ebinger et al., 2017), (2) it encompasses a broad range of spreading rates, crustal structures, and degrees of mantle melt , and (3) deformation has recently been detected at several calderas (Biggs et al., 2011) leading to new studies of their eruptive histories and geophysical properties Gislason et al., 2015;Lavayssiere et al., 2019;Martin-Jones et al., 2017;Samrock et al., 2015). Seismological constraints on crustal properties suggest the presence of vertically extensive intrusive complexes (Ebinger et al., 2017). Numerical models show that repeated intrusions into extensional settings can create mush zones in the upper and lower crust which remain partially molten for residence times of 10 5 -10 6 yr (Karakas & Dufek, 2015).
However, the processes and timescales by which upper mantle mafic and upper-crustal silicic magma reservoirs form and interact remain poorly quantified. The long (hundreds to thousands of years) recurrence intervals of eruptive activity and lack of high-resolution monitoring capabilities at individual volcanoes mean that geophysical estimates of the magma fluxes feeding silicic calderas are challenging.  Biggs et al., 2011) and historical eruptions (Tullu Moje, ∼1900; Kone, 1820, and Fentale, 1810, Wadge et al., 2016). The contours show the absolute shear velocity in the upper mantle, which indicates the presence of partial melt in focused upwellings along the Central Main Ethiopian Rift axis (Gallacher et al., 2016). (b) Map of gravity benchmarks at Corbetti and stack of Sentinel 1 interferograms from ascending path 160 (29 September 2015 to 17 October 2016). Each color cycle represents 2.8 cm in the satellite line-of-sight direction. (c) Vertical displacement at continuous GPS site CBM2 (syn. CURG Lloyd, Biggs, Birhanu, et al., 2018) relative to a stable reference, HAWA, located at the University of Hawassa showing a total uplift of 14.7 ± 2.7 cm from 31 March 2013 to 31 December 2016. The cGPS sites at Corbetti were installed in 2013.
Here we present joint ground deformation and time-lapse gravity data to elucidate the timescales, rate, and structure of magma supply beneath Corbetti caldera, an active silicic center in the Main Ethiopian Rift (MER) (Figure 1a). Although the spreading rate in the MER is relatively modest (Bendick et al., 2006) (6±2 mm yr −1 ), the rift is narrow, the melt supply is focused, and peralkaline rhyolites comprise 60-90% of erupted products (Abebe et al., 2007). Low seismic velocities in the upper mantle beneath the rift indicate that the presence of partial melt and absolute shear velocities at some points are lower than at some mid-ocean ridges (Figure 1a) (Gallacher et al., 2016). Combined with large population exposure and high vulnerability to environmental shocks, this makes understanding the silicic systems in the MER an important societal as well as scientific goal.

Caldera Deformation
Our ability to measure surface deformation at volcanoes remotely has dramatically increased the number of systems known to be deforming, and the diversity of deformation patterns observed. Periods of uplift lasting years to decades have been observed at silicic calderas around the world without resulting in immediate eruption (Phillipson et al., 2013), and indeed, deformation without eruption is statistically more common for volcanoes in caldera and rift environments (Biggs et al., 2014). Uplift alone cannot be uniquely attributed to fresh magma input, and the interpretation of associated monitoring signals (such as seismicity, gravity changes, and degassing) and the identification of eruption precursors remains a priority. Alongside improvements in monitoring, our understanding of the architecture of magmatic systems is rapidly evolving with new models on the incremental growth of crystalline mushes and images of trans-crustal magma systems (Cashman et al., 2017).

10.1029/2020GC008952
A new generation of numerical models aims to constrain the dynamics and evolution of sub-volcanic reservoirs by considering the thermomechanical feedbacks between magma and crust (Degruyter & Huber, 2014;Gottsmann et al., 2017;Karakas & Dufek, 2015). In a suitable thermal regime, magma influx can be accommodated inelastically without significantly increasing the reservoir pressure, yet still causing surface deformation (Degruyter & Huber, 2014;Head et al., 2019;Le Mével et al., 2016). Thermomechanical effects may reduce the amplitude and wavelength of surface deformation causing elastic models of surface deformation to overestimate source volume changes. The scientific debate on causes and consequences of volcano uplift has important implications for our understanding of the growth of large silicic magma bodies but until now has centered on a small number of well-studied systems in a limited range of tectonic environments (Hill, 2017).

Corbetti Caldera
Here we focus on Corbetti (Figure 1), a peralkaline caldera system formed at 182 ± 28 kyr, during a rift-wide ignimbrite flare-up, which would have dramatically altered the palaeoenvironment of the early hominids . Since ∼ 20 kyr two eruptive centers, named Chabi and Urji (syn. Wendo Koshe), have developed within Corbetti's collapse structure, dominantly composed of silicic deposits from pyroclastic falls and flows and aphyric obsidian lava flows. The youngest tephra layer, the Wendo Koshe Pumice, has been carbon-dated to ∼ 2.3 kyr b.p. (Rapprich et al., 2016), and terrestrial and lacustrine records suggest that an explosive eruption has occurred roughly once per 700-1,000 yr during the Holocene Martin-Jones et al., 2017). The Wendo Koshe Pumice is post-dated by four obsidian lavas, suggesting a return rate of 500-550 yr for effusive eruptions. Since 2009, Corbetti has been experiencing steady caldera-wide uplift at a maximum rate of 6.6±1.2 cm yr −1 , accompanied by microseismicity (Lloyd, Biggs, Birhanu, et al., 2018). Although over 500,000 people live within 20 km of Corbetti and geothermal power is being developed, all volcano geophysical studies have involved temporary scientific networks, and there remains no permanent, real-time volcano monitoring (Brown et al., 2015).

InSAR and GNSS
Sentinel-1 synthetic aperture radar (SAR) data from August 2015 to January 2017 were processed using the GAMMA software package (Werner et al., 2000) within the LiCSAR facility (González et al., 2015). While our data represent a subset of the 2007 to 2017 data considered in (Lloyd, Biggs, Birhanu, et al., 2018), we follow their processing approach. We select scenes with temporal baselines less than 48 days (20 ascending and 21 descending) to produce 66 interferograms. Final interferograms were resampled to 100 m. We produce a chain of ascending and descending interferograms, converted to rate, starting and ending at the closest interferograms to the dates of the gravity survey (ascending: 29 September 2015 to 17 October 2016, descending: 06 October 2015 to 30 September 2016) to use in the geodetic source modelling. We also processed one ascending and one descending Cosmo-SkyMed (CSK) interferogram (24 January 2014 to 19 January 2015 and 25 January 2014 to 11 December 2014, respectively) using the ISCE software package (Rosen et al., 2012). These dates were selected to match the period of the gravity studies. CSK interferograms were iteratively resampled and filtered to maximize coherence, with a final pixel size of 120 m following Lloyd, Biggs, Birhanu, et al. (2018) (Figure 1b).
We processed the coordinates of continuous GPS sites in the Corbetti area using GAMIT/GLOBK (Herring et al., 2010) in order to produce daily site position estimates with standard corrections for the ionosphere and modeled wet troposphere. To do so we used 15 International GNSS System (IGS) reference sites including ADIS, which is located approximately 200 km to the north of the study area. Daily solutions were generated in the International Terrestrial Reference Frame 2014 (ITRF14) (Altamimi et al., 2016). Figure 1c reports the position of the site CBM2 with respect to the reference site HAWA, located in the nearby city of Hawassa. For a broader comparison of GPS and InSAR data from Corbetti we refer the interested reader to (Lloyd, Biggs, Birhanu, et al., 2018) where data covering the period between 2009 and 2016 are reported. Here we focus on the period of microgravity measurements (2014-2016) and use a combination of GPS and InSAR (from ESA's Sentinel-1 constellation and ASI's CosmoSkymed) data to provide constraints on the displacement field.

Time-Lapse (Dynamic) Microgravity
Gravimetric data were collected in late September/early October of 2014, 2015, and 2016 using a Scintrex CG-5 gravimeter (serial no. 572) at three benchmarks (Table 1 and Figure 1b) and related to a reference  Lloyd, Biggs, Birhanu, et al., 2018), and resulting data are used for deformation modeling only (see Figure 7).
(CBAS) located outside the area of interest. This ensured that the effects of seasonality including those from inter-annual variations in gravity due to ground water changes were minimized. The number of benchmarks we were able to install and maintain over the observation period was limited by accessibility, vandalism of survey sites, and field crew safety. Corbetti is difficult to access and lies in an area which underwent political and ethnic violence during the observation period, and therefore, careful planning went into the selection of field and reference sites. A medical clinic was chosen to site the latter. To confirm that any displacement at CBAS is smaller than the displacement uncertainty, we used multiple interferograms to measure the relative LOS changes between the GPS location and locations further away from the volcano, where deformation would not be expected. Benchmark occupations followed a strict protocol whereby several repeat measurements were conducted as part of a survey loop, which started and ended at the reference. At least 5 cycles of 45 second-long readings at 6 Hz were performed at each benchmark, and the number of cycles was adjusted appropriate to the required data precision of <2 μGal. Individual benchmarks were visited between three and six times per survey resulting in data repeatability of between 1 and 7 μGal at individual benchmarks. Because of the difficult terrain and resultant shaking of the instrument during transport, survey loops were kept as short as possible and were generally completed within less than 4 hr. The raw gravity data were reduced for tidal effects of both Earth and ocean components using the Wahr-Dehant (Dehant et al., 1999) and GOT99.2 (Ray, 1999) latitude-dependent models, respectively.
At Corbetti, measured gravity changes result primarily from contributions from (i) vertical surface deformation, (ii) mass redistribution by the deformation of host rock by reservoir volume changes, (iii) mass changes within the reservoir, and (iv) density changes by compressing encasing rocks due to reservoir expansion.
Since the purpose of this study is to quantify contributions from (iii), gravity data were corrected for contributions from (i), (ii), and (iv) to obtain residual gravity changes. For (i) we used the theoretical free-air gravity gradient (FAG) = −308.6 μGal m −1 (whereby 1 μGal = 10 −8 m s −2 ) based on a measured FAG of −310 ±10 μGal m −1 at CBAS. At CBM2, we use the vertical component from the cGPS site. GPS benchmarks at CBM1 and CBM3 were vandalized, so we use interferograms from ascending and descending tracks to calculate the vertical uplift at the closest point in space and time following the approach outlined in . We quantify contributions (i)-(iv) numerically using finite-element analysis (see below) noting that contributions from (ii) and (iv) are dependent on reservoir shape (Trasatti & Bonafede, 2008;Walsh & Rice, 1979) and crustal properties (Currenti, 2014) and are expected to be significant in a mechanically heterogeneous and viscoelastic crust.

Analytical Modeling
To characterize the source of the deformation we test two analytical elastic deformation source models (a point and dislocation source) inverting the sub-sampled line-of-sight InSAR data and three-component cGPS site CBM2 for the source parameters using a Bayesian inversion approach described in detail in Lloyd, Biggs, Birhanu, et al. (2018) and references therein. We report the optimal value for each parameter from

Numerical Modeling
The protracted eruptive record and available geophysical, geological, and petrological observations imply the presence of a long-lived (>180 ka) magmatic system beneath Corbetti within a thermally primed and mechanically heterogeneous middle and upper crust. This has two major implications for the joint modeling of geodetic and gravity time series: (i) a hot upper crust facilitates ductile deformation with implications for the temporal displacement of subsurface density boundaries and (ii) mechanical heterogeneity by crustal density variations controls subsurface strain partitioning. Neither can be accounted for in traditional analytical geodetic models which demand an elastic, isotropic, and mechanically homogeneous crust. Therefore, we use finite-element analysis (FEA) to solve jointly and simultaneously for stresses, strains, and gravity changes in a thermomechanically heterogeneous crust. Employing the optimal location parameters for the elastic point-source model and the source paramater-space exploration for the period 2009-2016 presented in Lloyd, Biggs, Birhanu, et al. (2018) as starting parameters, the FEA incorporates a physically plausible crustal rheology based on petrological and thermal considerations consistent with a silicic system within a mature magmatic rift Iddon et al., 2018;Keranen et al., 2009) as well as geophysical constraints on the crustal structure (Lavayssiere et al., 2019).
We develop a suite of 2-D axisymmetric models using COMSOL Multiphysics (v5.3) with the free surface (1,650 m a.s.l.) at z = 0 km, horizontal extent of 40 km in both x and y directions, and 41.5 km in the vertical (z). This domain represents the crust in the southern MER, where the Mohorovicic discontinuity (Moho) is at ∼ −40 km b.s.l. (Keranen et al., 2009). The top surface is free to deform, the lateral surfaces have roller boundaries, and the bottom surface is fixed ( Figure 2b). See Table 2 for parameters, values, and symbols.
We convert the strain-based model of a point source to a stress-based model of a finite sphere and implement time-dependent mechanical properties using a Standard Linear Solid (SLS) parameterization (Head et al., 2019), consisting of a Maxwell element (elastic spring and dashpot) and an elastic spring in parallel where the instantaneous elastic response of the medium is governed by the shear modulus G(z), and the viscous response by a characteristic relaxation time 0 = ∕G(z). We consider the fractional components ( 1 and 0 ) of the shear moduli split equally across the two branches of the SLS, whereby 1 = 0 = 0.5. is parameterized by the Arrhenius approximation = A×e H RT , where A is a constant, H is the activation energy, R is the gas constant, and T is the absolute temperature.
Thermomechanical boundary conditions and domain parameterizations are informed by available geophysical and petrological constraints. We estimate mechanical parameters (Figure 2) using the 1-D p-wave (v p ) velocity model from Lavayssiere et al. (2019) to derive the crustal density ( c ) and dynamic Young's modulus E d using the relationships (Brocher, 2005) where Poisson's ratio = 0.23 is derived from the ratio between compressional and shear wave velocities presented in Lavayssiere et al. (2019). Shear modulus G and bulk modulus K for the crust are derived from E d ( Table 2).
It is important to note that E d is greater than the static Young's modulus E s by a factor of 1.3 to 3 for temperatures and pressures relevant for the slow upper-crustal deformation investigated here (Cheng & Johnston, 1981;Gudmundsson, 2012;Zhan et al., 2019). Here, we assume E d = E s with the implication that the deduced pressure changes and reservoir mass changes are maximum estimates. To specify the derived mechanical parameters E d and c in our model we use piece-wise cubic interpolations (Figure 2a). We fix the temperature at the top surface to 293 K and at the base of the model to 1273 K to match the background geothermal gradient in the area (Keranen et al., 2009). Ignoring radiogenic heat production in the crust, this parameterization results in the crustal shear viscosity distribution shown in Figure 2b.
The intrusion is implemented as a spherical source given inferences on source shape reported above. We explore source center depths within the range derived by the elastic model (z = −6.6 to −7.2 km) including an optimal depth of −6.9 km and source radii between 1.5 and 2.2 km. The dimensions of the source reflect the lateral extent and depth of the conductive anomaly imaged at a depth of 5 km and below (Gislason et al., 2015) as a proxy for the shallow-most part of a partially molten crystal mush zone. They also accounti for an upper limit of permissible annual pressure increments 10 MPa which are thus below common reservoir failure criteria (Gudmundsson, 2012). Source radii are chosen such that for any given source depth, the source does not extent significantly above the seismic brittle-ductile transition depth of z ∼ −6.5 km beneath the center of the deformation (Lavayssiere et al., 2019).
Analysis of incompatible trace elements from the MER suggests 80-96 vol.% fractionation of basalts to generate peralkaline rhyolitic melt and thus requires cumulate volumes of at least three to five times and perhaps 30-50 times greater than the erupted silicic volumes (Gleeson et al., 2017). Volume estimates of post-caldera silicic eruptions at Corbetti range between 12 and 35 km 3  indicating ≫100 km 3 of refractory material at depth. In our FEM, the shallow plumbing system is represented by a truncated cone with its geometrical and thermomechanical boundary conditions prescribed according to available geophysical data as follows. The top of this domain lies at a depth of 0.2 km below the surface and is 4 km wide, representing the hydrothermal system imaged by magnetotellurics and drilling (Gislason et al., 2015). The base of this domain has a radius of 3 km and lies at z = −9 km matching the depth of the 10 Ωm resistivity isosurface (Gislason et al., 2015). Thermally, the bottom of this domain represents a hot, but fully crystalline intrusive complex (1273 K). The shape of the resulting shallow storage zone-defined here as the region within the 1173 K isotherm and z < −10 km-is similar to that produced by thermal models for other magmatic systems (Paulatto et al., 2012), and the volume is ∼175 km 3 .
In our model the thermal conditions of the intrusive complex are implemented at sub-solidus conditions with the implication that viscoelastic effects may be more pronounced if the intrusive complex is at a higher temperature than modeled here.
A typical 2-D axisymmetric joint gravity and deformation model has ∼50,000 mesh elements solving for ∼2.5 M degrees of freedom (DOF). Pressure increments are represented by a load applied to the source boundaries. The source has mechanical properties consistent with a magma reservoir as specified in Table 2. We approximate the observed quasi-linear uplift at Corbetti since 2009 (Lloyd, Biggs, Birhanu, et al., 2018) and increase in residual gravity by quarterly ( t = 0.25 yr) increments of both source pressure and source density from t 0 = 2009 onwards. The model calculates optimal pressure and density change time series to fit the observed deformation and gravity changes between 2014 and 2016 by solving for starting values (at t = 0) of two step functions as shown in Figures 3a and 3b. Mean annual volume changes (derived by volume integration over the expanding source ) are calculated from predicted changes between t = 5.25 yr (October 2014) and t = 7.25 yr (October 2016).
We simultaneously calculate gravity changes from the density variations caused by the concurrent displacement field u and the redistribution of mass in the subsurface by solving Poisson's differential equation (Currenti, 2014) for the gravitational potential g, whereby ∇ 2 g = -4 GΔ (x,y,z). Dirichlet boundary conditions for the gravitational potential at infinity are set to zero, and three source terms quantify the different contributions of subsurface density changes Δ (x,y,z) to the gravity changes -g z via Δ (x,y,z) = -u · ∇ 0 + ic ∇ · u.
The first term on the right-hand side quantifies the contribution of the displacement of density boundaries in a heterogeneous medium (at the free surface, crustal density discontinuities, and the reservoir wall) to the observed gravity changes, while the second term results from the change in density by the intrusion of new mass. The third term is the density changes resulting from volume changes in compressible surrounding rocks. The term i is composed of two components: (i) the density change resulting from the input of new mass into the reservoir at constant volume V and the compression of resident magma and (ii) the density change accompanying the volume expansion V of the reservoir. The former component is dependent on the compressibility of the reservoir r = c +1/ i × i ∕ P where c is the compressibility of the country rock, while the latter depends on the density change induced by the inflating reservoir and the replacement of mass surrounding the reservoir. Residual gravity changes are derived from g r (x,y,z) = − g z -U z , where ( U z ) is the free-air effect.  (Lloyd, Biggs, Birhanu, et al., 2018) of ∼6.6 ± 1.2 cm yr −1 . Incremental mass addition with density increases as shown in (b) reproduce mean residual gravity changes of ∼9 μGal yr −1 observed between 2014 and 2016 (d). The data shown are for a spherical source with a radius of 2 km located at z = −6.9 km (5.3 km b.s.l.). The red lines (in c and d) mark the period of interest for the current study.

Ground Deformation and Gravity Changes
The current phase of protracted uplift at Corbetti began in 2009 and has been detected by several independent satellite systems (Biggs et al., 2011;Lloyd, Biggs, Birhanu, et al., 2018) and follows a spatially distinct earlier period of deformation associated with pressure changes within a shallow (<1 km) fault-bounded hydrothermal reservoir . Caldera-wide quasi-concentric uplift has now been occurring at a near-constant rate (albeit minor fluctuations) for more than 10 yr, and interpolating over data gaps between satellite systems gives a cumulative vertical uplift of >60 cm (Lloyd, Biggs, Birhanu, et al., 2018). A continuous GPS network installed in 2013 provided independent confirmation of the steady uplift. The site closest to the center of deformation, CBM3 (syn. C01G Lloyd, Biggs, Birhanu, et al., 2018), recorded an uplift rate of 6.5 ± 1.3 cm yr −1 , while the longest-running station, CBM2 (syn. CURG Lloyd, Biggs, Birhanu, et al., 2018), recorded a vertical rate of 3.9 ± 0.7 cm yr −1 (Figure 1c). Surface deformation and residual gravity changes are reported in Table 1, including associated uncertainties. Gravity readings at CBM3 were significantly noisier compared to the other benchmarks during the measurement cycles (periods of seconds to minutes) as well as during reoccupations within and between survey loops (periods of hours). This benchmark is located in a hydrothermally active area as shown by magnetotelluric measurements (Gislason et al., 2015), and we cannot rule out contamination of the signal over a broad range of frequencies  (Rymer & Williams-Jones, 2000) with respect to the free-air gradient (FAG) and the Bouguer-corrected FAG (BCFAG; shown for Bouguer slab densities between 2,500 and 2,800 kg m −3 ) (Mickus et al., 2007). Data from inside the caldera (Table 1) plot in sector S3 of the U z versus Δg space and are indicative of mass addition and density increase in the source. Sector S1 indicates both mass and density decrease, while sector S2 represents mass increase and density decrease. (b) Spatiotemporal variations in gravity changes across the survey network. The red star marks the center of the surface deformation. (c) Temporal evolution of gravity changes at survey benchmarks relative to the gravity reference CBAS. (d) Mean annual gravity changes as a function of horizontal distance from the center of the deformation source. The red line marks the optimal solution while the shaded area marks the range of solutions for subsurface mass and density changes for the joint deformation and gravity change model reported in Table 1. Values given in (b)-(d) are residual changes after free-air correction (see section 3). from hydrothermal processes that are difficult to capture using annual surveys (Gottsmann et al., 2007). The data uncertainty reported for this benchmark includes the noise observed during the surveys over periods of seconds to days.

Figure 4a shows observed gravity changes (Δg) and contemporaneous vertical displacements (Uz) between
October 2014 and October 2016. The gradient Δg/Uz is indicative of the mass and density changes at depth (Rymer & Williams-Jones, 2000). Data along the free-air gradient (FAG) result from surface uplift and imply a density decrease, while data following the Bouguer-corrected free-air gradient (BCFAG) result from additional displacements of density boundaries, implying a mass addition. Data in either of the three sectors of the U z versus Δg space marked S1, S2, and S3 in Figure 4 arise from different permutations of mass and density changes (S1: mass and density decrease; S2: mass increase and density decrease; S3: mass and density increase). Data from Corbetti are indicative of mass and (minor) density changes in the source.
After correcting for the effect of U z on the measured gravity, statistically significant residual changes of 12 ± 5 μGal and 18 ± 6 μGal are observed inside the caldera at sites CBM2 and CBM3, respectively (Table 1 and Figure 4b). The temporal evolution of these changes shows a constant increase through 2015 and 2016 (Figure 4c).
The observation of positive gravity residuals implies an addition of material with higher mass density than that of surrounding rocks. Previous Bouguer gravity and seismic observations suggest that the upper crust of the central MER is 2,400-2,470 kg m −3 above 5 km b.s.l. (Lavayssiere et al., 2019;Mickus et al., 2007) and ∼2,670 kg m −3 below 5 km b.s.l. (Mickus et al., 2007). Thus, the gravity field increase with time in addition to the caldera-wide concentric ground deformation precludes the hypothesis that the uplift is solely driven by low density material such as hydrothermal fluids, which have densities of 600-1,150 kg m −3 at super-critical to subcritical conditions (Zhang & Frantz, 1987).
From RMS values and Akaike information criteria of the tested source geometries the point-source solution provides a better fit, and this finding is in line with results reported by Lloyd, Biggs, Birhanu, et al. (2018), where model fits and residuals for a longer observation period between 2009 and 2016 are reported.
The inferred geodetic source depth range of 6-8 km is equivalent to a confining pressure of between ∼130 and ∼200 MPa using the available velocity model for Corbetti (Lavayssiere et al., 2019). Fractionation depths in the central and southern Main Ethiopian Rift are typically shallow (∼100-150 MPa) (Gleeson et al., 2017), suggesting that magma may have undergone some fractionation upon reaching the source zone. The depth range corresponds to the upper part of a low resistivity (<10 Ωm) body (Gislason et al., 2015), which for a peralkaline composition typical of Corbetti (73 wt % SiO 2 ; 4.9 wt % Na 2 O), would require >2% partial melt (Pommier & Le-Trong, 2011). Thus, from the combination of geodetic depth estimates based on analytical models and residual gravity changes alone, we propose that the source is likely magmatic.

Magma Density Estimates
To estimate the possible source densities at Corbetti, we consider two end-member (most mafic and most silicic) magma compositions erupted at Corbetti: alkali basalt and peralkaline rhyolite (Table 4 and Figure 5).
For the rhyolite, we use the composition of the Wendo Koshe Pumice at Corbetti . No major and trace elements compositions for Corbetti's basalts have yet been published, so we compare estimates from the 1.6 Ma Bofa basalts near Aluto volcano  and a set of post-caldera basaltic scoria cones at Kone (Iddon et al., 2018). We use standard values of partial molar, volumes, thermal expansion coefficients, and compressibilities of the oxide components (Spera, 2000) and   Fontijn et al. (2018), , and Iddon et al. (2018) and Standard Parameters From Spera (2000) Peralkaline Note. Volatile estimates (marked by * ) are from preliminary melt inclusion studies (Edmonds, pers. comm.). than 0.1 kg m −3 and it is not necessary to specify precise P-T conditions to obtain reliable melt density estimates. While hybrid magmas between residual mantle-derived basaltic melt and anatectic crustal melt may form, we ignore their contribution due to the comparably low expected volume fraction (<6%) (Karakas & Dufek, 2015).  Spera (2000). Volatile estimates are from preliminary melt inclusion studies (Edmonds, pers. comm.). Source data are provided in Table 4.
Melt inclusion studies from the Main Ethiopian Rift suggest ∼5 wt.% H 2 O in rhyolitic magmas at depths of ∼5 km and ∼1 wt.% H 2 O and 2,000 ppm CO 2 in basaltic magmas (Edmonds, pers. comm.). Including the dissolved volatile components gives wet melt densities of 2,084 ± 60 kg m −3 for the Wendo Koshe Pumice, 2,623 ± 120 kg m −3 for the Kone Scoria Cones, and 2,595 ± 110 kg m −3 for the Bofa Basalt. Magma is likely to contain crystals and/or exsolved volatiles in addition to melt, which will increase or decrease the density, respectively. We estimate the density of a fully crystalline alkali gabbro cumulate  to be ∼3,000-3,100 kg m −3 and ∼2,700 kg m −3 for a peralkaline granite. Given that magma viscosity increases dramatically at crystal contents of 50% (Costa, 2005), we take this as an upper limit for mobile magmas. This gives a range of plausible magma densities of 2,000-2,450 kg m −3 for peralkaline rhyolites and 2,500-2,860 kg m −3 for alkali basalts. The presence of a free gas phase will lower these values significantly. The density values compare to crustal densities derived from seismic velocity data as shown in Figure 2.

Numerical Deformation Source Models
Due to the similarity in the footprints of the geodetic and microgravity signals (Figure 4b), we hypothesize a single source behind the observed surface uplift and gravity changes at Corbetti. We select a sphere as a geologically plausible approximation of Corbetti's upper-crustal magmatic plumbing system based on (1)   Note. We use latitude, longitude, depth, and volume change values from the best fit elastic geodetic inversion models as starting parameters. Derived pressure changes and radii reflect values matching an upper limit of source pressure changes of 10 MPa yr −1 and a source size such that it does not extend >100 m above the seismic brittle-ductile transition depth (Lavayssiere et al., 2019) at z ∼ −6.5 km. The mean mass flux is derived from V r × i , the lower and upper bounds are derived from the products between lower bound volumes and upper bound density changes, and vice versa, to account for the trade-off between source volume and density changes. and (2) geophysical imaging of kilometer-thick intrusive complexes beneath silicic centers within the MER and wider EARS and more specifically, magnetotelluric images of the sub-volcanic architecture at Corbetti (Gislason et al., 2015).
Best fit uplift velocities and gravity changes in the center of the deformation anomaly as a function of time from the numerical model are displayed in Figures 3c and 3d. For the period 2014-2016 the model predicts an average uplift rate of 6.6 cm yr −1 matching the observations (Lloyd, Biggs, Birhanu, et al., 2018) in combination with an average gravity change of 9 μGal yr −1 . Figure 6. Modeled contributions to observed annual gravity residual data from Corbetti (2014Corbetti ( -2016 given an inflating spherical source by mass intrusion. The contributions to the residual gravity changes include (i) the free-air effect given a gradient = −308.6 μGal m −1 , (ii) the displacement of density boundaries due to source expansion in mechanically heterogeneous host rock, (iii) the source density change, and (iv) the density change caused by the compression of host rock. The example shown is for a spherical source with a radius of 2.0 km, a source center at z = −6.9 km, and undergoing a density increase of 7 kg m −3 yr −1 .
The numerical models fit the deformation and gravity patterns for source radii of 1.9 ± 0.2 km and source depths of between 7.0 ± 0.2 km (Table 5 and Figures 6 and 7). Annual source pressure changes are predicted in the range between 1 and 10 MPa. Modeled source volume changes (0.88 ± 0.04 × 10 7 m 3 yr −1 ) for the period 2014/2016 are lower than those from the analytical models (Table 3) which assume a mechanically uniform and elastic crust. We find that the combination of a mechanically heterogeneous compliant upper crust and ductile deformation results in a source volume change of ∼1/3 lower than the value proposed by Lloyd, Biggs, Birhanu, et al. (2018), whereby crustal compliance contributes to about half of this reduction in volume change. The derived annual source density changes are 7 ± 4 kg m −3 , which for background densities in the range ∼2,400-2,670 kg m −3 (Lavayssiere et al., 2019;Mickus et al., 2007) is consistent with either dry peralkaline rhyolite magma or mobile alkali basaltic magmas ( Figure 5 and Table 4). Modeled source volume and source density changes are inversely related with a resultant source mass change of between 1.11 and 1.95 × 10 11 kg yr −1 .
We derive a reservoir compressibility r in the range between 0.1 and 1.0 GPa −1 for magma densities between 2,400 and 2,670 kg m −3 , annual pressure changes between 4 and 10 MPa, and density changes between 3 and 10 kg m −3 for the best fit models ( Figure 8). The compressibility values are indicative of a magma reservoir containing compressible fluids in the presence of a free gas phase. Estimates of magma compressibility range from 0.02 GPa −1 for unsaturated or degassed magmas to >10 GPa −1 for porous magmas (Huppert & Woods, 2002;Rivalta & Segall, 2008).

Discussion and Implications
Data interpretation from joint geodetic and microgravity studies is often challenging because the inferred source density change is sensitive to assumptions on the source geometry and the subsurface mechanics. Both aspects are difficult to resolve with limited spatial data coverage and the use of traditional homogenous, isotropic, and elastic half-space models. To capture the wider deformation field and the complex subsurface mechanics and resultant density variations at the restless Corbetti caldera, this study makes use of two areas of technological advances: (1) satellite images to constrain source geometry, which even downsampled, provides over a thousand measurements of displacement along two independent line-of-sight vectors, and (2) finite-element analysis which considers a thermomechanically heterogeneous crust consistent with long-term magmatic activity and the effect of crustal mechanical heterogeneity and an inelastic rheology on the mass redistribution by deforming host rock.
The combination of high-resolution satellite data with high-precision time-lapse gravity data and numerical models provides a fresh perspective on a long-lived debate regarding the mechanisms of caldera unrest in line with the current understanding of the architecture of long-lived magmatic systems.
The residual gravity increase at Corbetti unequivocally demonstrates the influx of magma beneath the caldera, whereby a significant portion of the mass flux is being accommodated inelastically in the upper crust and within a compressible reservoir. The derived reservoir compressibility of between 0.1 and 1 GPa −1 Figure 8. Parameter space for the best-fit joint deformation and gravity change models. Panel (a) shows the reservoir compressibility as a function of reservoir density, annual reservoir density change, and annual reservoir pressurization. While density has only a minor influence on compressibility, higher compressibility values imply higher annual density changes and lower annual pressurizations. Reservoir compressibility in the best fit models varies by 1 order of magnitude (see color bar). Panel (b) gives the dependence of reservoir compressibility on the deduced annual magma flux rate for given annual reservoir density changes. The gray area marks best-fit reservoir density change values of between 3 and 10 kg m −3 yr −1 . The best-fit annual magma flux rate ranges between 1.1 and 1.95 × 10 11 kg yr −1 for a 1 order of magnitude change in reservoir compressibility.
suggests the presence of a free gas phase, and we propose that up to 1,000 ppm CO 2 (and negligible quantities of H 2 O) could be exsolved from alkali basalts at the source depth (Papale et al., 2006). The implication of a compressible reservoir and low Vp/Vs ratios suggest that at least some of this exsolved gas remains trapped within the magma reservoir (Lavayssiere et al., 2019).
Our approach clearly demonstrates the importance of considering inelastic behavior and magma compressibility when modeling ground deformation and gravity change time series at long-lived magmatic systems and, by doing so, provides physically plausible solutions. Employing thermomechanical constraints, we produce a numerical model that explains the data by modest annual density increases and pressure increments below common reservoir failure criteria (Gudmundsson, 2012). Our preferred source model behind the protracted uplift at Corbetti is the intrusion of unfractionated alkali basalt into Figure 9. Sketch of the upper-crustal plumbing system beneath Corbetti caldera derived from the combination of the geodetic and gravimetric surveys and data modeling in conjunction with other geophysical observations highlighting the seismic brittle-ductile transition zone (SBDTZ after Lavayssiere et al., 2019) at z ∼ −6.5 km and the 10 Ωm resistivity isosurface (deep conductor after Gislason et al., 2015). In this interpretation, the current magma intrusion is within a partially molten crystal mush zone. A shallower peralkaline rhyolite plumbing system which fed the explosive and effusive post-collapse eruptions at Corbetti is expected to reside at the top of the intrusive complex from which it fractionates. A hydrothermal system at 0.2-1 km below the surface represents the main geothermal reservoir beneath the caldera (Gislason et al., 2015) but is not shown for clarity. The underpinning thermal model includes a background geotherm after Keranen et al. (2009) with a temperature of 1273 K at the Mohorovicic discontinuity (Moho) while a thermal aureole represents an intrusive magmatic complex within which the source is embedded.
an upper-crustal reservoir containing a free gas phase ( Figure 9). Nonetheless, insufficient geophysical constraints exist to fully constrain a numerical model for Corbetti's, or indeed any, magmatic system, and models should be refined with new constraints from subsurface imaging or physical properties as they become available.
The mass flux indicates that intrusive complexes beneath silicic caldera systems can undergo recharge at a much higher rate than previously considered for mature continental rifts. The deduced flux is up to 2 orders of magnitude higher than the mean mass eruption rate at Corbetti over the past 70 ka  and supports petrological models of pantellerite formation and geophysical constraints of underpinning mafic magmatism: (i) 90-96% fractional crystallization of alkali basalt is required to produce peralkaline rhyolite melt (Gleeson et al., 2017) which contributes 60-90% of eruptive products in the MER (Abebe et al., 2007) and (ii) a middle to lower-crustal high seismic velocity zone along the rift axis interpreted to be fully crystalline high density intrusives (Ebinger et al., 2017;Mackenzie et al., 2005).
Extrapolating the 2014-2016 reservoir volume change to the 10 yr period of steady uplift gives a cumulative source volume change of 0.09 ± 0.02 km 3 for Corbetti ( Figure 10). Although the magnitude and duration of the current uplift at Corbetti is unique in the East African Rift during the available survey period, several other peralkaline caldera systems are known to be deforming (Biggs et al., 2011). Many of the signals are consistent with multi-year periods of magma input and subsequent readjustment of the hydrothermal system , suggesting that the combined magma flux into silicic systems along the rift may be spatially and temporally variable, but consistently high.
In Afar and Iceland, where the spreading rate is considerably faster, the magma volumes derived for basaltic mega-dyke intrusions (>2 km 3 ) (Wright et al., 2012) are at least one of magnitude greater than the cumulative volume change observed at Corbetti (∼ 0.1 km 3 ). However, the Corbetti value matches or exceeds the magma volumes associated with most recent dyke intrusions in regions of slow to moderate extension  (Xu et al., 2015), and Saudi Arabia in 2009 (0.13 km 3 ) (Pallister et al., 2010), and shows no sign of declining (Lloyd, Biggs, Birhanu, et al., 2018). In the Main Ethiopian Rift, dyke intrusions have been postulated based on far-field geodetic transients (Bendick et al., 2006), but none have yet been directly detected suggesting that they are small, slow, deep, or infrequent (Biggs et al., 2011). Thus, our observations suggest that the time-averaged flux associated with upper-crustal intrusions into silicic reservoirs is comparable to and likely higher than that associated with short-lived basaltic dyke intrusions.
We suggest that substantial volumes of mafic magmas are periodically supplied and fractionated beneath long-lived silicic rift caldera systems to form and sustain partially molten ("mushy") and compressible upper-crustal intrusive complexes with new intrusions accommodated largely by ductile deformation. The intrusion responsible for the ongoing uplift at Corbetti represents a significant component of the upper-crustal magma budget in the MER and is a modern analog that illustrates the rapid processes by which magmatism and tectonics modify the lithosphere to accommodate strain in protracted magma-assisted mature continental rifting. Higher-resolution three-dimensional geophysical images are required to determine whether upper-crustal intrusive complexes are isolated beneath individual magmatic centers of the MER or extend laterally along the rift axis. Silicic eruptions are triggered by rapid destabilization and remobilization of partially molten intrusive complexes (Allan et al., 2012;Burgisser & Bergantz, 2011), and this may have implications for hazard and risk assessment at Corbetti and the wider MER.

Author Contributions
J. B. and J. G. conceived the project, J. B., J. G., R. L., and Y. B. conducted the field work and data collection, Y. B. processed the GPS data, R. L. processed the InSAR data and performed analytical geodetic modeling. J. G. processed the gravity data and developed the numerical models. J. B. and J. G. wrote the initial manuscript draft, and all authors contributed to the writing of the final manuscript.