What Triggers Caldera Ring-Fault Subsidence at Ambrym Volcano? Insights From the 2015 Dike Intrusion and Eruption

where Abstract Surface deformation accompanying dike intrusions is dominated by uplift and horizontal motion directly related to the intrusions. In some cases, it includes subsidence due to associated magma reservoir deflation. When reservoir deflation is large enough, it can form, or reactivate preexisting, caldera ring-faults. Ring-fault reactivation, however, is rarely observed during moderate-sized eruptions. On February 21, 2015 at Ambrym volcano in Vanuatu, a basaltic dike intrusion produced more than 1 m of coeruptive uplift, as measured by InSAR, synthetic aperture radar correlation, and Multiple Aperture Interferometry. Here, we show that an average of ∼ 40 cm of slip occurred on a normal caldera ring-fault during this moderate-sized (VEI < 3) event, which intruded a volume of ∼ 24 × 10 6 m 3 and erupted ∼ 9.3 × 10 6 m 3 of lava (DRE). Using the 3D Mixed Boundary Element Method, we explore the stress change imposed by the opening dike and the depressurizing reservoir on a passive, frictionless fault. Normal fault slip is promoted when stress is transferred from a depressurizing reservoir beneath one of Ambrym's main craters. After estimating magma compressibility, we provide an upper bound on the critical fraction ( f = 7%) of magma extracted from the reservoir to trigger fault slip. We infer that broad basaltic calderas may form in part by hundreds of subsidence episodes no greater than a few meters, as a result of magma extraction from the reservoir during moderate-sized dike intrusions. its own weight. Caldera faults, or cracks surrounding the reservoir which extend from the reservoir to Earth's surface, form as the lid collapses. Ambrym volcano (Vanuatu) has a 12-km wide caldera, and researchers propose it formed during an explosive eruption 2,000 years ago. However, in 2018, Ambrym's caldera sunk along caldera faults during a nonexplosive eruption. This observation questions whether an explosive eruption was necessary to form Ambrym's caldera in the first place. Furthermore, in February 2015, an eruption 10 times smaller than in 2018 also caused the ground to sink along caldera faults. Utilizing ground motion data obtained from satellite radar systems to model magma reservoir outflow and fault displacement, we conclude that, in 2015, the ground sank along caldera faults. This sinking is explained by the removal of as little as <7% of the stored magma from the reservoir. We therefore propose that Ambrym's wide caldera may have formed as a result of many frequently occurring medium-sized eruptions. This challenges the thought that wide calderas mainly form as a result of large eruptions.

faults break up the collapsing piston (e.g., Bárðarbunga [Ágústsdóttir et al., 2019], Kīlauea [K. R. Anderson et al., 2019;Neal et al., 2018;Segall et al., 2020], and Piton de la Fournaise [Peltier et al., 2009]). In contrast, volcanoes that host shallow and broad magmatic chambers tend to collapse in a coherent fashion, and the central portion of the collapsing piston stays intact (Roche et al., 2000). Walker (1984) hypothesized that silicic calderas may form incrementally by downsagging, which can be accompanied by caldera ring-faulting. The fraction of material removed from a magma reservoir needed to activate caldera ring-faults does not depend strongly on eruptive style. Incremental caldera ring-fault activation may also occur during frequent eruptions at broad, shallow basaltic calderas. Notably, at Sierra Negra, which hosts the largest caldera in the Galápagos (9 km maximum diameter), there is no geological evidence of a catastrophic caldera collapse (e.g., lithic breccias and ignimbrite deposits [Druitt & Bacon, 1986]) in the literature, despite voluminous historical eruptions (Munro & Rowland, 1996;Reynolds et al., 1995).
It is challenging to confirm that incremental subsidence is an important driver of caldera development at broad, shallow basaltic systems. Any significant subsidence would only accumulate over time scales of hundreds of years. Little geological trace would be left at the surface, especially if lateral intrusions, which may arrest at depth, were the dominant mechanism of magma withdrawal from a central magmatic plumbing system. Increased space-geodetic monitoring of volcanoes worldwide improves the chances of measuring caldera subsidence and ring-fault activation during moderate-sized eruptions (VEI < 3) (Pinel et al., 2014;Shreve et al., 2019).
Caldera ring-faults can be activated when the reservoir pressure drops beneath lithostatic (i.e., is underpressurized) during magma extraction (Druitt & Sparks, 1984). At basaltic caldera-rift systems, this extraction occurs primarily through lateral dike intrusions (Sigmundsson, 2019). However, in some cases of dike propagation, regardless of whether or not ring-faults have been activated, surface displacements associated with a depressurizing reservoir are not observed geodetically. This may be due to a masking of the subsidence signal by large displacements related to the dike intrusion (Grandin et al., 2009). In addition, other factors such as host rock or magma compressibility can reduce the amplitude of the geodetic signal (Rivalta & Segall, 2008). Accounting for mechanical source interactions allows additional constraints to help identify a depressurizing source.
In February 2015, a dike intrusion fed an eruption at Ambrym volcano and produced more than 1 m of asymmetrical uplift to the SW of the fissure, as measured by synthetic aperture radar (SAR) data. To understand how the dike intrusion led to caldera subsidence and ring-fault reactivation, we combine 3D Mixed Boundary Element Method (BEM) calculations with a neighborhood inversion algorithm to determine the deformation sources contributing to the surface displacement. We then calculate the static stress change on the ring-fault to investigate the possible contribution of the dike opening and reservoir deflation on caldera ring-fault reactivation, as well as determine bounds on the pressure and volume change of the magma reservoir needed to trigger caldera ring-faulting.

Ambrym Volcano, Vanuatu
Ambrym is a basaltic volcanic island located in the central portion of the New Hebrides subduction zone in Vanuatu (Figure 1a). The island hosts two rift zones oriented N105°S, which radiate bilaterally from the island's central 12-km wide caldera (McCall et al., 1970) (Figure 1b). The caldera has been previously interpreted as resulting from the collapse of a giant tuff cone during a sequence of explosive phreatomagmatic eruptions (Robin et al., 1993). This view was later challenged by Cronin and Németh (2005), based on fieldwork reporting the presence of altered mafic deposits around the island, but a lack of dacites. In December 2018, the first geodetically monitored rift zone intrusion occurred, with >0.4 km 3 of magma traveling more than 20 km into the SE rift zone, resulting in a submarine eruption that emitted basaltic pumice onto the nearby shoreline (Shreve et al., 2019). Coeval with the 2018 rift zone intrusion, the caldera floor subsided by more than 2 m. The caldera ring-faults also reactivated, most notably in the northern half of the caldera. A depressurizing deformation source was modeled at between 1 and 5 km beneath the summit, depending on how many sources were included (Hamling et al., 2019;Shreve et al., 2019). This event confirms that the caldera ring-faults of Ambrym are part of an active fault system, not solely relict structures resulting from a major collapse 2,000 years ago (Robin et al., 1993). Due to these observations of meter-scale caldera subsidence, in conjunction with historical documentation of reoccurring rift zone intrusions, recent studies (Hamling et al., 2019;Shreve et al., 2019) have invoked the hypothesis that Ambrym's caldera developed as a result of hundreds of similar rift zone intrusions, which simultaneously caused the gradual, incremental subsidence of the caldera floor. The rift zone intrusion in 2018 provided an example of caldera ring-fault reactivation at Ambrym. However, caldera faults may reactivate during smaller events. More frequently occurring moderate-sized eruptions, such as an intracaldera fissure eruption that occurred in 2015, provide constraints on the magma volume extraction necessary for ring-fault reactivation.
Ambrym's volcanic activity over the past decades includes lava lakes and Strombolian explosions within the nested pit craters located in the two main volcanic cones, Marum and Benbow, near the western caldera rim (see Figure 1c) (Németh & Cronin, 2008). In addition to lava lake activity, occasional intracaldera fissure eruptions have occurred, most notably in 1986 and 1988-89 (Eissen et al., 1991) (Figure 1c). The former took place in the eastern portion of the caldera, and erupted lava with a slightly elevated SiO 2 content (60 wt% SiO 2 ), compared to eruptive products from the main cones (∼50.5 wt% SiO 2 ). This suggests various degrees of melt differentiation within the magmatic system (Eissen et al., 1991;Robin et al., 1993).   Allard et al. (2015). (c) Enlargement of Ambrym's caldera, including the volcanic cones of Marum and Benbow, caldera rim, historical volcanic vents, and dates corresponding to historical lava flows. Lighter shades of gray indicate older volcanic deposits, and the lava flow from February 2015 is shown in orange. (d) Difference between a preeruptive TanDEM-X DEM and a posteruptive Pléiades DEM ( Figure S1). Newly opened fissures are denoted by thin dotted lines. (e) The cumulative seismic moment and daily SO 2 mass (as measured by the Ozone Mapping and Profiler Suite, OMPS) emitted throughout the duration of the eruption, which is marked by the orange dotted lines. The timing and focal mechanism of the 6.4 M w earthquake on February 19, 2015 is shown in green. Bounds on the SO 2 mass estimate are shown in light blue and are discussed in the Supporting Information Text S1.

February 2015 Eruption
On February 20-21, 2015, another intracaldera fissure eruption occurred, with a main fissure located ∼3 km SE of Marum ( Figure 1d) Hamling & Kilgour, 2020). In addition to the main fissure to the SE of Marum, a secondary fissure opened closer to Marum, ∼200 m south of Niri Mbwelesu Taten, oriented N133°S, which fed a 900 m-long lava flow (see Figure 1d). This event was preceded by a 6.4 M w earthquake on February 19, 2015, at 13:18 UTC with a hypocenter 30 km southeast of Ambyrm's craters and a depth of ∼5 km (Figures 1a, 16.50°S, 168.28°E, according to the location from the Oceania Regional Seismic NETwork, event ID ird2015dmnz). The earthquake had a reverse focal mechanism and moment tensor solution that included only 61% double couple component (Figure 1a). Thermal anomalies were detected by both MODVOLC (R. Wright, 2016;R. Wright et al., 2004) and the Middle Infrared Observation of Volcanic Activity system . According to , the eruption lasted for ∼44 hours, initiating sometime between 20 February 14:30 UTC and 21 February 2:40 UTC and ending, at the latest, on 22 February 11:10 UTC ) (see Supporting Information Text S1 for confirmation of this onset time using the Hybrid Single-Particle Lagrangian Integrated Trajectory model [Crawford et al., 2016;Stein et al., 2015] and information on the SO 2 plume height derived from Infrared Atmospheric Sounding Interferometer [IASI] observations [Clarisse et al., 2014]).  convert the middle infrared radiance into a time averaged lava discharge rate, estimating a total of 4.8 (±2.4) × 10 6 m 3 of lava was emitted from the main fissure oriented ∼N135°S, traveling ∼3 km.
To directly constrain the lava flow volume, we calculate a digital elevation model (DEM) using MicMac software (Rupnik et al., 2018) and posteruption Pléiades optical satellite images with no cloud cover over the lava flow. We estimate an updated lava flow volume of ∼12.4 (±0.08) × 10 6 m 3 , or ΔV erupted = ∼9.3 (±0.08) × 10 6 m 3 DRE, with an average flow thickness of ∼5 m (Figure 1d, Supporting Information Figure S1, Text S2, and Table S1). Uncertainties were calculated using the methods of Bagnardi et al. (2016) and Favalli et al. (2010). This value is ∼2 times larger than the value estimated by .
We use the Ozone Mapping and Profiler Suite (OMPS) observations to conservatively estimate ∼40 kt of sulfur dioxide (SO 2 ) emitted during the eruption (Figure 1e, Supporting Information Text S3). For this SO 2 mass estimation, we used the OMPS TRM product, which assumes a center of mass altitude of ∼7.5 km . This is in agreement with the retrievals of the SO 2 plume height from IASI observations (see Text S3). An SO 2 mass estimate of ∼40 kt is consistent with degassing of a lava volume corresponding to the ∼9.3 × 10 6 m 3 lava flow DRE (see Supporting Information Text S3). After February 25, SO 2 degassing returned to passive background levels (∼7 kt/day, according to Carn et al., 2017). These estimates of erupted lava volumes and emitted gas will be compared with volume changes at depth, as constrained by geodesy. This comparison will provide constraints on the mass balance between the erupted and intruded material, and the overall size of the magma reservoir at depth.

Geodetic Data
In this study, we exploit multiple SAR data sets to measure the coeruptive ground displacement (Table S2). A joint analysis of differential InSAR and Multiple Aperture Interferometry from ALOS-2, and pixel offset tracking from COSMO-SkyMed (CSK), allows us to decompose the 3D coeruptive displacement field. We use the 3D displacements to interpret the deformation signal, determining the possible location and nature of the active deformation sources.

Differential InSAR
The Japan Aerospace Exploration Agency's L-band SAR satellite ALOS-2 acquired images before and after the eruption, in both ascending (stripmap) and descending (stripmap and ScanSAR) orbit geometries. With a wavelength of 24 cm, interferograms produced using data acquired from L-band SAR satellites maintain coherence in vegetated regions, such as tropical volcanic islands. We process interferograms from ascending stripmap mode (SM3) images spanning January 24 to March 21, 2015 with the Interferometric SAR scientific computing environment (ISCE) (Rosen et al., 2012). Filtering and unwrapping are performed with NSBAS modules (Doin et al., 2011;Grandin et al., 2012;Rosen et al., 2004). After geocoding, the final pixel posting is ∼14 × 14 m. Supporting Information Text S4 describes in detail the processing steps.
Before the eruption, the ALOS-2 descending archive consists of only ScanSAR (interferometric wide-swath mode, WD1) acquisitions. Coherent ScanSAR-to-ScanSAR interferograms can only be calculated with images acquired after February 8, 2015, due to inadequate (less than 50%) burst synchronization on the ground for acquisitions before this date (Lindsey et al., 2015;Natsuaki et al., 2016). This was due to an issue with ALOS-2's navigation system, which was fixed on February 8, 2015 (Lindsey et al., 2015). As a result, standard ScanSAR-to-ScanSAR processing is replaced by ScanSAR-to-stripmap processing for descending images spanning January 17, 2015 to March 14, 2015 (Kobayashi et al., 2015;Natsuaki et al., 2016;Ortiz et al., 2007). This is implemented with the GSISAR software developed by the Geospatial Information Authority of Japan (Kobayashi et al., 2015;Tobita, 2003). A 30-m mesh Advanced Spaceborne Thermal Emission and Reflection Radiometer Global Digital Elevation Model was used to remove the topographic fringes (Tachikawa et al., 2011). The postinterferogram formation steps (multilooking, filtering, unwrapping, and geocoding) are performed with NSBAS modules, as described in the Supporting Information Text S4. The descending interferogram is multilooked one time in range and eight times in azimuth, resulting in a final pixel posting of ∼30 × 20 m after geocoding.

Multiple Aperture Interferometry
In addition to differential InSAR, which measures satellite line-of-sight (LOS) displacements, we also process a multiple aperture interferogram (MAI) (Bechor & Zebker, 2006;Liang & Fielding, 2017), in order to derive the along-track, or azimuth, displacement using ascending SM3 images spanning January 24 to April 4, 2015. These measurements are then used in the 3D displacement decomposition. Azimuth common-band filtering with a normalized squint of 0.66 generates subaperture single-look-complex pairs which are used to calculate forward and backward looking interferograms with a multilook factor of 8 and 16 in range and azimuth, respectively. The difference between these interferograms provides the along-track displacement, albeit with a lower accuracy than differential InSAR. The postprocessing steps include filtering (Goldstein & Werner, 1998), unwrapping, and geocoding with a pixel spacing of ∼85 × 70 m.

Pixel Offset Tracking
Subpixel offset tracking measures surface displacement by finding the cross-correlation peak of two SAR image amplitudes (Michel et al., 1999). Pixel offsets are less accurate than InSAR (Michel et al., 1999). However, they are complementary as they measure both LOS and azimuth displacement, as well as provide measurements in areas where phase is not coherent or surface displacements are large.
For CSK descending acquisitions spanning 13-25 February 2015, we calculate pixel offsets in both range and azimuth. We also calculate pixel offsets from the ALOS-2 ascending pair spanning January 24 to March 21. All pixel offset calculations are run using ISCE with a window size of 64, a skip width of 32, and a search width of 20. Images are then geocoded, and postprocessing includes scaling by pixel size (3.6 m/pixel for ALOS-2 azimuth, 1.1 m/pixel for CSK range, and 2.1 m/pixel for CSK azimuth), as well as referencing to the median of a 15 × 15 pixel box to the northwest of the caldera (near 168.05°E, 16.24°S). We also remove unreasonable displacement values, filter, mask based on the signal-to-noise ratio (SNR), and mask densely vegetated areas outside the caldera (Table S3).

3D Decomposition
Following the method of T. J. , we invert for the 3D coeruptive displacement field using the InSAR LOS measurements, pixel offsets, and MAI measurements. This results in an overdetermined system of equations, and we perform a least squares inversion to solve for the horizontal (NS and EW) and vertical components of displacement (see Figure 2). Each data set was weighted equally in the inversion, which assumes similar accuracy for all measurements. InSAR has a higher accuracy than the other data sets, but the equal weighting ensures that contributions from the pixel offsets and MAI, such as displacement along the satellite azimuth, are not underprioritized in the 3D decomposition (Grandin et al., 2018). Notable characteristics of the 3D displacement field allow us to hypothesize which deformation sources were active during this event.

Coeruptive Displacement Field Description
The ascending and descending interferograms measure, respectively, up to 1.5 and 1 m of LOS shortening (ground motion toward the satellite). The 3D decomposition reveals that this discrepancy is due to significant horizontal motion in the displacement field. The cross sections in Figure 2 display the ratio of vertical to horizontal motion. In addition to horizontal motion, we emphasize three notable characteristics of the coeruptive 3D displacement field: 1. There exists a N-S asymmetry and discontinuity across the eruptive fissures in all components (see Fig-ure 2, Profile A-A′). The maximum LOS shortening is located just to the south of the main fissure, and this signal is most likely associated with the dike intrusion feeding the eruption. The ratio of subsidence to the north of the fissures compared to uplift to the south implies a dike dipping toward the SW. 2. There is also an E-W asymmetry in the western portion of the caldera, across a region that follows the caldera rim, and is incoherent and discontinuous in the interferograms (see Figure 2, Profile B-B′). This asymmetry is also observed in all components. The ratio of horizontal to vertical displacement indicates that this signal may be attributed to normal slip on a portion of the caldera ring-fault. As measured by CSK pixel offsets, this displacement occurred before 6:00 UTC on February 21, during the onset of the eruption ( Figure S3). 3. The 3D displacement field also estimates >20 cm of subsidence in the coherent region to the NW of the dike, beneath Marum crater ( Figure 2). The ascending interferogram shows ∼24 cm (two fringes) of LOS lengthening in the northern portion of the caldera, following the northeastern caldera rim. However, this vegetated portion of the caldera is incoherent in the descending interferogram, limiting the ability to exploit two independent measurements to confirm whether these fringes are attributed to ground displacement or atmospheric effects. This signal, confined to within the caldera, is reminiscent of deformation that occurred in the month following the 2018 rift zone intrusion. The 2018-2019 deformation was modeled by a shallow (4-5 km depth) depressurizing sill (Shreve et al., 2019).

Inversion Scheme
The displacement field's complexity hints at contributions from multiple sources. Geodetic modeling of volcanic deformation often capitalizes on computationally efficient analytical solutions to describe deformation due to pressure sources (Mogi, 1958), as well as shear or tensile dislocations (Okada, 1985). When multiple sources are involved, the displacement fields produced by each of these sources in an elastic, homogeneous medium are added linearly. Any mechanical source interaction, or stress transfer, between the sources is implicitly neglected.
Given the hypothesis that multiple sources are at play during the 2015 Ambrym eruption, and the geological complexity of caldera systems (preexisting weaknesses and faults due to caldera collapse [Acocella, 2007]), we choose to proceed with the 3D Mixed BEM (see Section 4.2 for more details) (Cayol & Cornet, 1997, 1998. This numerical method calculates the surface displacement due to pressure changes within volume sources (or "massive boundaries") or within fractures with complex geometries (such as curvature, varying dip, etc.), as well as the mechanical source interaction between all modeled sources. By combining the 3D Mixed BEM forward models with a nonlinear inversion, we are able to estimate the sources' geometries, while simultaneously inverting for a uniform stress change on the structures (i.e., pressure or shear stress changes).

Data Subsampling and Covariance Matrix
To compare modeled and measured surface displacements, measurements must first be subsampled to include a reasonable number of datapoints in the inversion, while still retaining enough information to robustly estimate the relevant model parameters. We include and subsample both the ascending and descending InSAR measurements in the inversion. Neither the SAR pixel offsets, nor the MAI displacements measurements, are included due to their low SNR. However, we project modeled surface displacements into the LOS for all available data sets a posteriori ( Figure S4). The data are subsampled using an adaptive quadtree decomposition algorithm (Jónsson et al., 2002;Walstead, 1999), as described in the Supporting Information Text S5, Table S4, and Figure S5.
Displacement values at subsampled points may be spatially correlated due to phase contributions from atmospheric, ionospheric, or other noise sources (Sudhaus & Jónsson, 2009). To account for this spatially correlated noise in the interferometric phase, the correlation distance and data variance are used to populate a data covariance matrix, according to , where  2 d is the variance, a c is the correlation length, and ‖k 1 , k 2 ‖ is the Euclidian distance between two subsampled pixels k 1 = (x 1 , y 1 ), and k 2 = (x 2 , y 2 ) (Tarantola, 1987). Following Fukushima et al. (2005)  3 10 d m 2 and a c = 330 m. The data covariance matrix weights the data in the inversion, taking into account the correlated noise between two pixels.

3D Mixed Boundary Element Method
As previously mentioned, we estimate the geometry and stress change of multiple pressure sources with a 3D Mixed BEM numerical approach (Cayol & Cornet, 1997, 1998. In the 3D Mixed BEM forward calculation, boundary conditions are the traction-free topography and uniform pressure changes on either surfaces (tensile cracks or shear fractures) or massive boundaries (reservoirs or the ground surface). Triangular elements are used to discretize all boundaries. The surface topography is derived from a 12-m resolution Tan-DEM-X DEM ( Figure S6a) (B. Wessel, 2016). Stress changes can either be imposed on the meshed triangular elements of the fractures or massive boundaries, or result from stress transfer from nearby sources, or both.
The 3D Mixed BEM combines the displacement discontinuity and direct displacement methods (Crouch, 1976;Lachat & Watson, 1976;Sokolnikoff, 1956). A description of these methods can be found in the Supporting Information Text S6. Nine parameters control the fracture geometry ( Figure S7) (Fukushima et al., 2010). The fractures are connected to the topography mesh by user-defined surface fissures. In this case, two sets of fissures are organized in an "en echelon" geometry. While their location has a degree of uncertainty, we use high-resolution optical satellite imagery from the Pléiades constellation to map the fissures that opened during the eruption (see Figure 1d). For the caldera ring-fault, the surface trace is determined using the SAR pixel offsets, in conjunction with identifying the incoherent regions in the InSAR measurements.
Displacements and stresses are discretized and interpolated on the boundaries. Thus, the more refined the mesh, the more precise the solutions. There is a trade-off between accuracy and computational efficiency of the forward model, controlled by the number of meshed elements used to represent the boundaries (see Figure S6b for a plot of model run duration vs. accuracy). In order to have a reasonable accuracy versus efficiency balance, the mesh is finer close to the echelons, where the displacement gradient is largest, and coarser in the far-field, where little-to-no displacement is measured ( Figure S6a). A shear modulus of 2 GPa and a Poisson's ratio of 0.25 are chosen, consistent with in situ measurements and depth of the inverted structures (Cayol & Cornet, 1998) and laboratory measurements (Heap et al., 2019). The influence of the chosen shear modulus on the estimated pressure change and reservoir volume will be discussed further in Section 6.2.3.

Nonlinear Inversion: Neighborhood Algorithm
The stress change and geometry of the fractures or reservoirs are inverted using the neighborhood algorithm-a simple, derivative-free search method (Fukushima et al., 2005;Sambridge, 1999a). The neighborhood algorithm utilizes nearest neighbor regions (i.e., Voronoi cells) to optimize the parameter space search and obtain an ensemble of models that best fits the data. An initial exploration of the parameter space picks N 1 random combinations of the estimated model parameters. To ensure a sufficient exploration of the initial parameter space, N 1 increases exponentially with the number of inverted parameters (Sambridge, 1998;Tridon et al., 2016). Each subsequent iteration picks N 2 combinations of parameters, and the N 2 computations can be performed in parallel (Fukushima et al., 2005;Smittarello et al., 2019). We take N 2 = 48 to optimize this parallelization, which was performed on a cluster using three nodes, with 16 CPUs per node. A particu-lar iteration searches within the N 2 multidimensional Voronoi cells with the lowest misfit from the previous iteration.

The misfit is defined as
where u o is the data, u m is the modeled displacements, and C d is the data covariance matrix, as defined in Section 4.1. This iterative procedure continues until a threshold criterion, σ NLST , is reached, based on the standard deviation of the misfit values from the last NLST forward calculations. We set NLST = N 2 . To calculate posterior probability density functions (PPDs), marginal PPDs, mean model parameters, and model uncertainties, we follow the Bayesian inference framework (Tarantola, 1987) described in Sambridge (1999b) and implemented by Fukushima et al. (2005) (Supporting Information Text S7).
The inversion procedure inverts for a phase offset but does not invert for ramp parameters in the InSAR data that may be due to residual orbital errors (Zebker et al., 1994). Therefore, we mask the displacement signal, fit a ramp to the nondeforming areas, and subtract the ramp from the ascending interferogram using General Mapping Tools (P. Wessel et al., 2019). For the descending interferogram, nearly the entire coherent region is deforming, so we do not estimate and remove a ramp.

Final Inversion Strategy
Using the inversion scheme outlined previously, we run three inversions of increasing complexity to test the necessity of adding multiple deformation sources. The number of nonlinearly inverted free parameters and number of forward models for each inversion are listed in Table 1. In order to justify increasing complexity (i.e., adding a new source), we calculate the Akaike Information Criterion (AIC) for each inversion, which takes into account the number of inverted parameters (i.e., given a particular misfit, the model with less inverted parameters is preferred). The AIC is defined as AIC = 2k + χ 2 + log|C d | + N × log(2π), where k is the number of inverted parameters, χ 2 is the misfit, |C d | is the determinant of the data covariance matrix, and N is the number of data points (Akaike, 1974). In each inversion, we use the same subsampled data sets and data covariance matrix from Section 4.1. Therefore, to compare the AIC of inversions, we can simply compare AIC = 2k + χ 2 .

Dike (Inversion 1)
According to (1) in Section 3.3, the N-S asymmetry across the eruptive fissure implies a dike dipping toward the SW. The first inversion estimates 10 parameters nonlinearly-the nine dike geometry parameters and the pressure change. Table 2 lists all inverted, fixed, and postprocessed parameters. The best fit model and data, synthetics, and residuals are shown in Figure 3. After the inversion is complete, the best fit model dike is meshed with smaller elements. The final synthetic surface displacements are scaled by a constant to account for discrepancies between computations performed with the coarse and fine meshes (Fukushima et al., 2005). This scalar is ∼1 for all inversions (see Table S6, indicating no significant difference in the surface displacements calculated using the coarse and fine meshes). The mean model and 95% uncertainty intervals, as well as the marginal PPDs, are listed in the Supporting Information for all inversions (Table S5 and Figures S9-S11). The best fit model dike has a scaled pressure change of ΔP dike = 2.2 MPa and the total volume change is ΔV dike = 24 × 10 6 m 3 . Volume change is calculated by integrating the opening on each element, which itself is determined by the source geometry and pressure change (Figures 3a and 4a; Table S6).

Dike and Fault (Inversion 2)
The residuals from Inversion 1 reflect more than 20 cm of unexplained displacement in the western portion of the caldera, corresponding to the signal discussed in (2)  Note. Each inversion ended when the standard deviation criterion (<0.5) was reached. both a dike and a frictionless fault, whose top edge is predefined and intersects the topography mesh. We invert for three additional parameters-the dip, bottom edge depth, and a uniform shear stress change on each fault element, resulting in normal fault slip (i.e., stress change imposed parallel to the orientation of the fault elements). To reduce the number of nonlinearly inverted parameters and ensure computational feasibility, we do not invert the angle of the bottom edge of the fault, the twist, shear, or any curvature (see Figure S7 for a description of these parameters). This simplification is chosen because analog and numerical models have found caldera ring-faults to be vertical at depth (Beauducel et al., 2004;Gudmundsson et al., 2016;Liu et al., 2019). We also only include the portion of the fault that we hypothesize was reactivat-SHREVE ET AL.

10.1029/2020JB020277
10 of 26 Note. Green cells indicate that the parameter is inverted, red cells indicate that the parameter is fixed, orange cells are the parameters that are calculated after the inversion, and gray squares indicate that the source is not included in that inversion. See Figure S7 for an explanation of the listed parameters and Table S5 for the explored intervals.   ed in 2015 (see Section 6.1). Because the fault and dike orientations are perpendicular, pressurization of the dike results in closing (fracture wall interpenetration) on some sections of the fault plane ( Figure S8). We impose a nonnegativity constraint to avoid fracture wall interpenetration. This nonnegativity constraint is enforced in the 3D Mixed BEM through Lagrange multipliers, as described in Cayol et al. (2014) (from hereon we call this constraint "preventing interpenetration"). As a result, computation time of a single forward model increases from 30 s to 11 min for the source geometry shown in Figure 4b.
The dike opening and fault shear displacements of the best fit model are shown in Figure 4b. The dike's scaled pressure change is 2.7 MPa, the dike's total volume change is ΔV dike = 23 × 10 6 m 3 , the fault's external shear stress change is −0.21 MPa, and the fault's average slip is 0.39 m. The AIC decreases by more than 17% in Inversion 2 compared to Inversion 1 (Table 1 and Figure 3b).

Dike, Fault, and Reservoir (Inversion 3)
There remain long-wavelength features in the Inversion 2 residuals, postulated to result from reservoir depressurization, with a similar spatial footprint to the caldera floor subsidence measured after the Ambrym 2018 rift zone intrusion (see [3] in Section 3.3). Therefore, a final inversion includes a thin (semiminor axis of 2.5 m), horizontal, oblate spheroid structure, whose location (beneath Marum crater, at 168.12°E, 16.25°S) and depth (4.1 km b.s.l.) are fixed according to the postintrusion deflation source estimated in Shreve et al. (2019). Inverting for these parameters would result in an unreasonably long inversion. As shown by Sambridge (1998), the number of initial computations N 1 needed to adequately explore the parameter space increases exponentially with the number of inverted parameters. Test inversions also show that the depth cannot be adequately constrained, and we therefore rely on independent data sets (e.g., Allard et al., 2015;Hamling & Kilgour, 2020;Hamling et al., 2019;Legrand et al., 2005;Shreve et al., 2019, for reservoir depth) and postprocessing (see Section 6.2.2) to obtain reasonable values. Again, we invert for normal and shear stress changes on the dike and fault as nonlinear parameters. We approximate the reservoir as a thin oblate spheroid as opposed to a horizontal fracture because the nonnegativity constraint restricts closing on all fractures in the inversion. Instead of closing of a horizontal fracture, we impose a uniform negative pressure change boundary condition on the spheroid. For a given ΔV res , reservoir volume V and change in pressure ΔP res are inversely proportional (Amoruso & Crescentini, 2009), resulting in the high ΔP res obtained by the inversion given the thin reservoir. After running 3D Mixed BEM forward models for spheroids with various thicknesses, we find that the difference in both the surface displacement and the source volume change is within reasonable error bounds for all aspect ratios (see Section 6.2.2 and Figure S12).
We invert for the same parameters as Inversion 2, as well as the reservoir's horizontal radius (assuming it is axisymmetric) and pressure change. The source opening and shear displacements of the best fit model of Inversion 3 are shown in Figure 4c, the marginal PPDs are shown in Figure S11, and the residuals are shown in Figure 3c. Pressure change (scaled) of the best fit model dike is 2.1 MPa, the total volume change of the dike is ΔV dike = 24 × 10 6 m 3 , the total volume change of the reservoir is ΔV res = −15 × 10 6 m 3 , and the average closing on the reservoir is −0.9 m. The fault's external shear stress change is −0.13 MPa and the average slip on the fault is 0.44 m. Fault slip reaches ∼1 m at depth, close to the reservoir. We note that the AIC lowers 39% between Inversions 1 and 3 (see Table 1).

Inversion Comparisons
The final estimated dike geometries diverge primarily at depth (Figure 4), due at least in part to the lack of sensitivity of the surface displacement to pressure changes from deep sources (Du et al., 1992). The mean fault dip estimated in both Inversion 2 and 3 is within a few degrees, however, the dip tends toward a shallower value. We decide to constrain the minimum dip to 65° based on numerical and analog modeling, as well as geological observations, of steep, inward-dipping normal caldera ring-faults for calderas with low roof aspect ratios (R a , ratio of chamber depth to width of caldera surface expression) , and references within). According to geodetic modeling from the 2018 eruption ( , but whether or not this case is relevant for Ambrym is beyond the scope of this study. We acknowledge that the similarities between the fault dip for Inversions 2 and 3 may be artificially imposed by the parameter limits. The total volume change in the dike is 23-24 × 10 6 m 3 for the three inversions. All estimates and uncertainties for pressure changes, external shear stress changes, total volume changes, and maximum and average slip values can be found in Tables S5 and S6. 6. Discussion

Reservoir Depressurization-Not Diking-Induces Faulting
In Section 5.3, we note that Inversion 3 has the lowest AIC, indicating that it is the most likely source configuration. We conclude that contributions from three sources-a dike intrusion, caldera ring-fault, and reservoir-best fit the data, even after the AIC takes into consideration the additional inverted parameters.
We will now discuss whether the configuration of three deformation sources is more mechanically consistent than only a dike and ring-fault, providing a supplementary argument, in addition to the AIC, for the presence of a reservoir.
In Inversion 2, we invert for external shear stress change on the caldera ring-fault in order to fit the surface displacements. We consider this an "external" shear stress change because an external stress perturbation must generate the shear stress change, which occurred either prior to the eruption (i.e., the fault was prestressed), during the event itself, or some combination of these events. Static stress transfer from the dike opening during the eruption is the most obvious source of perturbation.  calculated by τ CSC = τ + μ′σ, where τ is the shear stress change, σ is the normal stress change (positive is unclamping), and μ′ is the effective friction coefficient, where μ′ = μ (1 − B). B is Skemptons coefficient, which relates pore pressure to confining stress, and we set μ′ = 0.4, approximated from laboratory values (Nostro et al., 1998). A positive τ CSC indicates that the stress perturbation promotes fault failure. Because we have no information regarding the initial stress state, Coulomb stress change calculations can only be used to estimate if the fault has been brought closer to failure, not if the yield strength has been reached.
We calculate the Coulomb stress change in 3D on normal dip-slip fault planes, coinciding with the fault elements estimated in the inversions. The source geometries are fixed based on the results of Inversions 2 and 3. In the former, the normal and shear stress changes due to a pressurizing dike are computed using the Mixed BEM. For each fault mesh element i, the strike, dip, normal and shear stress change are used to calculate τ CSC,i . The first row in Figure 5 shows the Coulomb stress change on each element of the fault from Inversion 2. We find that the dike opening inhibits normal dip-slip at the fault's southern end, where surface displacements are dominated by compression perpendicular to the dike plane (Rubin & Pollard, 1988). Normal dip-slip is also inhibited at depth. The Coulomb stress change calculation from Inversion 2 has a slight increase in the northern portion of the fault, in front of the dike where the stress state is dominated by extension. We then calculate the Coulomb stress change from Inversion 3, when the stress perturbation results from a combination of a pressurizing dike and a depressurizing reservoir. There is a positive Coulomb stress change at the fault's north end, both at the surface and at depth, which promotes normal fault failure (see Figure 5).
Numerical modeling has shown that the stress field of a caldera is complex, influenced by the combination of regional tectonic stresses, edifice loading, unloading due to caldera formation, and other sources of local stresses (Buck et al., 2006;Corbi et al., 2015;Pinel & Jaupart, 2000). Given the unknown initial stress state of the caldera, we cannot discount the possibility that the fault was brought closer to failure due to stress transfer from only the dike intrusion, given the slight CSC increase. However, as shown in Figure 5, the Coulomb stress change calculations indicate that fault reactivation at depth is mainly promoted by the reservoir's depressurization as it fed the dike intrusion and fissure eruption, rather than by stress transfer from the dike intrusion itself. In fact, in Inversion 2, stress transfer from the dike intrusion causes a maximum of 0.21 m of reverse fault slip at the surface, inconsistent with measured surface displacements. In addition, if we calculate the Coulomb stress change due to only a depressurizing reservoir on a ring-fault encircling the entire caldera (assuming a vertical dip, given the lack of knowledge on ring-fault geometry beyond the western portion of the caldera), normal dip-slip along the entire western rim of the caldera is promoted (see Figure S13b). When the same calculation is performed with only a pressurizing dike, normal dip-slip is promoted in regions of the western rim where we do not observe surface displacement associated with faulting (see Figure S13c).
Nevertheless, stress transfer from the depressurizing reservoir alone is not sufficient to explain the magnitude of fault slip corresponding to the measured surface displacements. The stress transfer reaches a maximum value of 0.5 MPa on the fault ( Figure 5). In Inversion 3, we accordingly invert for an external shear stress change on the fault, representing either prestress from a previous event or stress transfer from a perturbation during the eruption. This external shear stress change may result from the accumulation of previous reservoir depressurization events, and the stress transfer from the dike and reservoir during the 2015 eruption brought the fault to its failure threshold. Between Inversions 2 and 3, the mean external shear stress change on the fault required to explain the magnitude of slip decreases by 38% (from −0.21 to −0.13 MPa, see Figure 5). In other words, in Inversion 3, less external "prestress" is required to produce the observed slip on the caldera fault, because the depressurizing reservoir provides a subsequent fraction of this stress perturbation. This indicates that the depressurizing reservoir induces fault shear displacements which produce deformation that are congruent with the observations.
The ring-fault's initial stress state plays an important role in estimating the critical fraction of magma withdrawal. In Inversion 3, the external shear prestress is −0.13 MPa, which is approximately equivalent to the shear stress transferred by the reservoir. We hypothesize that each fissure eruption and reservoir depressurization event of this size may induce ∼−0.13 MPa of shear stress on the fault. A previous withdrawal of magma from the reservoir, with about the same volume as the 2015 event, may have stressed the ring-fault, but not to the point of failure. We note that a nonzero friction coefficient would increase the magnitude of required external prestress, further biasing our estimates. Critical volume fractions will therefore be (at least) doubled, to account for the preexisting external stress introduced in the inversion.
However, there may be other physical processes that depressurized the reservoir and brought the fault closer to failure, accounting for the −0.13 MPa of preexisting external stress on the fault. For example, Ambrym's persistently high rates of passive degassing over the past decade may have caused reservoir depressurization, if there is no magma replenishment from a deeper source (Girona et al., 2014). Therefore, the critical volume fractions estimated below are upper bounds, because smaller volumes of magma extraction may cause ring-fault reactivation if the ring-fault was primed by other reservoir processes unrelated to magma extraction.

The Influence of Caldera Roof Aspect Ratio
Caldera ring-faults are formed during caldera collapse. In basaltic systems, calderas develop due to lateral magma propagation which drains a central plumbing system, activating caldera ring-faults and resulting in gradual collapse of the caldera (Sigmundsson, 2019). The plumbing system's geometry influences the critical fraction of magma needed to trigger collapse (defined as where Δq f is the magma extraction volume before the onset of collapse and V is the reservoir volume). Analytical, analog, and numerical models conclude that thresholds are lower for shallower, broader calderas (low roof aspect ratios, where h is reservoir depth and r is the reservoir radius) (Geshi et al., 2014;Holohan et al., 2011;Roche et al., 2000). Studies of caldera collapse events at Kīlauea, Piton de la Fournaise, Fernandina, Miyakejima, and Bárðarbunga conclude that f crit calculated from observations is smaller than estimated from analog modeling (K. R. Anderson et al., 2019, and references within). Ambrym has an especially low roof aspect ratio (<0.4) (Hamling et al., 2019;Shreve et al., 2019), estimated by comparing the depth of the main depressurizing deformation source in 2018 and the width of the caldera surface expression. The 2015 caldera ring-faulting event at Ambrym, although not a full-scale caldera collapse event, allows us to draw bounds on the critical volume fraction needed to reactivate preexisting caldera ring-faults at a broad, shallow basaltic caldera.

Pressure Change and Minimum Reservoir Volume
Using the 3D Mixed BEM and independent observations, we adjust the reservoir dimensions to establish bounds on ΔP res and V, after fixing the shear modulus μ (here we have μ = 2 GPa, see Section 6.2.3 for a discussion of uncertainties). These bounds will allow us to estimate magma compressibility and constrain f, which we here define as the critical fraction of magma needed to reactivate (not form) caldera ring-faults.
The volume V of an ellipsoidal reservoir depends on the semiaxes a, b, and c. Because trade-off exists between ΔP res and V (the full relationship for various reservoir geometries can be found in Amoruso & Crescentini, 2009). Here, c represents the "height" or "thickness" of the reservoir, which we fix in the inversion because of this trade-off. We explore the range of reasonable V by setting the semiminor axis c of the reservoir estimated in Inversion 3 to values that range between 2.5 and 1,300 m, adjusting ΔP res to obtain a ΔV res ≈ −15 × 10 6 m 3 (as estimated in Inversion 3). The residuals between the forward 3D Mixed BEM calculations for various aspect ratios and the synthetic displacement field from Inversion 3 are shown in Figure S12. Values are within 5% error of the synthetic displacement field from Inversion 3, and we conclude that we cannot further constrain the geometry of the reservoir tapped during the 2015 eruption. However, we note that a sill-like reservoir (a ≫ c) best explains previous episodes of caldera floor subsidence, such as the subsidence measured in 2015-2017 (Shreve, 2020) and after Ambrym's 2018 rift zone intrusion (Shreve et al., 2019).
We can put a lower bound on the reservoir's negative pressure change, because we do not observe a fullscale caldera collapse creating new ring-faults. A large underpressure ΔP − could cause the formation of new ring-faults and a full-scale caldera collapse. Using geometrical arguments of  (assuming vertical faults) and a maximum host rock shear strength τ max of 100 MPa (Schön, 1996;Touloukian et al., 1989), the largest underpressure for the estimated R a of Ambrym (∼0.4) would be −ΔP − = −4R a τ max = −160 MPa. We will take −ΔP max ≈ −175 MPa, assuming that the reservoir overpressure before the eruption was at most 15 MPa, given reasonable rock tensile strengths (K. Anderson & Segall, 2013;Delgado et al., 2019;Sparks, 1997). Given the estimated reservoir dimensions of a ≈ 1,340 m from Inversion 3, to obtain a ΔP − < 175 MPa, c can range from 5 to 1,300 m. This results in ΔP res ranging between −4.1 and −78 MPa and the corresponding V ranging from 0.02 to 9.5 km 3 . Given this range, ΔV res varies from −13 to −15 × 10 6 m 3 .
Although the estimated V ranges over multiple orders of magnitude, we note that Allard et al. (2015) estimate Ambrym's minimum reservoir volume to be 0.5 km 3 . Allard et al. (2015) obtain this volume by multiplying the magma residence time (estimated from radionuclide fluxes and activity ratios) by the estimated magma influx into the shallow magma reservoir, which sustains the measured SO 2 flux. To find this minimum, Allard et al. (2015) assume that if the degassed melt is not erupted at the surface or intruded into the crust, then it has to be recycled into a deeper portion of the magmatic system. On the other hand, to find a maximum estimate of V, we use passive degassing estimates from Carn et al. (2017) (7 kt day −1 ) and neglect magma recycling into a deeper magmatic system. In this case, the entire magmatic system would contribute to degassing over a period of ∼10 years, the time period of the Carn et al. (2017) study. Neglecting effects of magma compressibility, we estimate that 6 km 3 of magma was degassed, given ∼25,550 kt of SO 2 emissions over 10 years (see Supporting Information Text S3 for an example calculation). This provides an independent estimate of V.
Therefore, the independent bounds obtained from geodetic modeling, radioactive equilibria, and SO 2 flux span a range of overlapping values <10 km 3 (0.02-9.5, 0.5, and 6 km 3 , respectively). In addition, estimates for the reservoir volume change during the 2018 rift zone intrusion lie between 0.3 and 0.7 km 3 (Hamling et al., 2019;Shreve et al., 2019). If we assume that the reservoir volume change was 0.7 km 3 in 2018 and that the entire reservoir did not drain, we conclude that the shallowest compartment of Ambrym's reservoir (i.e., the portion that contributes to feeding lava lake activity and degassing) has a minimum size of ∼1 km 3 , yet upper bounds on the reservoir size cannot be robustly constrained in this study.

Magma Compressibility
By comparing the ratio r v of reservoir volume change to the volume of intruded material, we can estimate to the first-order the magma compressibility (Rivalta, 2010;Rivalta & Segall, 2008). Magma compressibility is one factor that contributes to the discrepancy between estimates of reservoir volume change and the amount of magma extracted from the reservoir (Johnson et al., 2000). Rivalta and Segall (2008) and Rivalta (2010) derived the generic formula for r v based on chamber compressibility β c and magma compressibility β m , such that where ΔV s and ΔV c are the volume change in the sink (dike + erupted volume) and reservoir, respectively. β c expresses how the host rock responds elastically to pressure change,  and an aspect ratio c a . We calculate chamber compressibility using analytical expressions from Amoruso and Crescentini (2009), finite element calculations (K. Anderson & Segall, 2011), and the 3D Mixed BEM. When a = b ≫ c (a penny-shaped sill), using the analytical expressions derived in Amoruso and Crescentini (2009), the following equation can be used to calculate chamber compressibility: Although Equation 2 fails for shallow penny-shaped sills (K. Anderson & Segall, 2011), the "medium" depth (4-5 km b.s.l.) of Ambrym's reservoir allows us to reasonably approximate β c using this expression for lower values of the aspect ratio c a . For the upper bound, chamber compressibility can be approximated by a sphere within 25% error when aspect ratios c a range between ∼0.6 and 10 (K. Anderson & Segall, 2011), such that Using Equations 2 and 3, it follows that β c ranges from 3.8 × 10 −10 to 22 × 10 −10 Pa −1 . The lower bounds are based on a spheroidal geometry, and the upper bounds are based on an assumption of a penny-shaped sill with an aspect ratio of ∼0.1 (or c ≈ 130 and V ≈ 1 km 3 , as concluded in Section 6.2.2). We confirm these analytical approximations by calculating β c using the 3D Mixed BEM. β c in this case ranges from 3.7 × 10 −10 to 18 × 10 −10 Pa −1 . The differences between the analytical and the computational values are consistent with K. Anderson and Segall (2011), and we will proceed with the discussion using the latter range of values.
From Equation 1, given ΔV dike and ΔV res estimated in Inversion 3 and ΔV erupted from Section 2.2, then . It follows that magma compressibility ranges from 4.5 × 10 −10 to 22 × 10 −10 Pa −1 . Although spanning almost 1 order of magnitude, we can conclude that this value is higher than that for degassed basalts (Spera, 2000). This estimate is consistent with recent studies of the Ambrym 2015 eruption, which hypothesize that Ambrym's shallow reservoir was pressurized due to bubble nucleation and growth (Hamling & Kilgour, 2020). The existence of gas bubbles within the chamber may explain the relatively high magma compressibility.
Estimates of magma compressibility at Kīlauea during the 2018 lava lake drainage and rift zone eruption ranged from 2 × 10 −10 to 15 × 10 −10 Pa −1 , also higher than values estimated assuming gas-poor basalt (K. R. Anderson et al., 2019;Segall et al., 2020). K. R. Anderson et al. (2019) and Segall et al. (2020) conclude that this may indicate the presence of bubbles in the reservoir, even though the value of β m relies strongly on prior estimates of the shear modulus μ and . Similarly, at Erta Ale, another volcano which hosts a lava lake, no significant cumulative postintrusion subsidence was measured after the 2017 intrusion, indicating either highly compressible magma or magma fed from a deep source (Moore et al., 2019). The observations listed above regarding high compressibility magma at Ambrym, Kīlauea, and Erta Ale need to be investigated in more detail before further conclusions can be drawn regarding magma compressibility at volcanoes hosting active lava lakes.
The magma compressibility estimate at Ambrym implies that the reservoir feeding the eruption was not fully degassed. Magma may instead be degassed at shallow levels through persistent lava lake activity. This is consistent with Allard et al. (2015), who concludes that Ambrym's degassing occurs in a "closed system" (Edmonds, 2008). "Closed system" degassing at Ambrym is consistent with the mass balance between the degassing during the eruption and the volume of erupted material. If the melt degasses at very shallow levels, we could assume that there would be no significant degassing through the dike. When comparing the degassed SO 2 mass, conservatively 40 kt, which is equivalent to a degassing of 9.5 × 10 6 m 3 of lava, we indeed notice this is the same order of magnitude as the volume of erupted lava (see Section 2.2). Therefore, we conclude that none of the degassing during this eruption was related to degassing from deeper within the magmatic system, whether from the dike intrusion or the reservoir. The same conclusion was drawn during the December 2018 eruption (Shreve et al., 2019). This balance between the degassed sulfur mass and erupted lava volume has already been documented at other effusive volcanic eruptions (Barnie et al., 2016), yet this is not necessarily the case in explosive eruptions, as emphasized by the "excess" sulfur problem (Kilbride et al., 2016;Wallace, 2001).
In the above discussion, we assume a shear modulus of 2 GPa, a value that is consistent with recent laboratory studies estimating upscaled values of elastic moduli for volcanic rocks (Heap et al., 2019). However, the host rock elastic moduli of a specific region are affected by factors such as porosity, temperature, the presence of microcracks, or confining pressure (Heap et al., 2019). These factors, among others, may cause the shear modulus to vary over more than an order of magnitude. Therefore, we acknowledge that, with no data to constrain the crustal properties at Ambrym, the shear modulus may range from 0.2 to 20 GPa. As a consequence, the chamber compressibility (which is inversely proportional to the shear modulus) may be as low as 3.7 × 10 −11 Pa −1 or as high as 1.8 × 10 −8 Pa −1 . Magma compressibility would then vary over a similar range. Comparisons with other magma compressibility estimates should be applied with caution, until studies at Ambrym can further constrain the shear modulus (such as with experimental rock strength estimates using local samples [Heap et al., 2019], or with seismic wave velocities [Albino et al., 2018;Ellis et al., 2007;Grandin et al., 2010;Masterlark et al., 2012;Wauthier et al., 2015]).
An alternative explanation for the discrepancy between the ΔV s and ΔV c is magma recharge of the plumbing system. However, even in the months following the 2018 rift zone intrusion, no posteruptive uplift related to magmatic recharge was measured, despite SAR acquisitions every 6 days (Shreve et al., 2019). In addition, the caldera floor subsided in the months following the 2015 eruption, as measured with both ALOS-2 and Sentinel-1 (Figures S14 and S15). Therefore, at this point, no evidence for magma replenishment can be found.

Critical Volume Fractions for Ring-Fault Reactivation
Given the estimated β c , β m , V, ΔP, and ΔV, we determine the critical volume fraction needed to trigger caldera ring-faulting. Only the western portion of the caldera ring-fault is reactivated, because the reservoir drained during the 2015 event is located beneath the active vents and lava lakes, and this is thus where the largest Coulomb stress change occurred on the ring-fault (See Figure S13b). In the 2018 rift zone intrusion and reservoir drainage event, the northern portion of the caldera ring-fault was reactivated, possibly due to withdrawal of a larger quantity of magma from a more central storage area within Ambrym's magmatic system (Hamling et al., 2019;Shreve et al., 2019).
After doubling the volume change estimates from Inversion 3, as explained in Section 6.1, we assume a volume of reservoir change −ΔV res,f = 30 × 10 6 m 3 , a volume of extracted material −Δq = r v ΔV res,f = 67 × 10 6 m 3 , and a minimum reservoir volume of 1,000 × 10 6 m 3 . Because the reservoir deflation and caldera ring-fault reactivation occurred before the end of the eruption ( Figure S3), the estimated Δq at the time of ring-fault reactivation is at most the sum of the magma volume intruded in the dike and the volume erupted at the surface. The upper bound on the critical fraction   Δ crit V V V is thus at most 3%, and, after taking into account the compressibility of the system, the critical fraction   Δq f V is at most ∼7%. These values are consistent with the fraction of extracted material needed to initiate fracture formation at the surface f exp , as found in analog models. Geyer et al. (2006) estimate that f exp may be as low as 2% for calderas with the same roof aspect ratio as Ambrym (r ∼ 0.4). The estimates for V crit and f, however, are upper bounds, considering Δq and ΔV may be overestimated, and the total volume of stored magma at Ambrym may be significantly larger than the estimated reservoir volumes given in this study (see Section 6.2.2).

Implications for Basaltic Caldera Formation
The 2015 eruption at Ambrym is an example of caldera ring-fault reactivation during a moderate-sized eruption. It is possible that localized caldera subsidence may occur during moderate-sized eruptions, which occur on decadal time scales at Ambrym. Caldera-wide subsidence events, on the other hand, occur during rift zone intrusions, because they typically involve a greater volume of magma withdrawal (Shreve et al., 2019). According to Eissen et al. (1991), more than 10 rift zone eruptions have occurred at Ambrym in the past 120 years. Assuming fault slip similar to that measured during the 2018 rift zone intrusion (∼0.4 m), this would result in ∼4 m of fault slip in a little over 100 years. Assuming the time frame estimated by Robin et al. (1993) since caldera formation (∼2,000 years ago) and stability of Ambrym's activity, we calculate almost 80 m of fault slip related to reservoir drainage due to rift zone intrusions and eruptions. The largest subsidence events at Ambrym, which occur during rift zone intrusions, contribute the most to caldera development. As observed in 2018, a more laterally extensive and slightly deeper magma lens may be tapped during these events, compared to the moderate-sized eruption in 2015. Better constraining the scaling relations between (1) intrusion size and frequency and (2) intrusion size and fault slip magnitude, would allow for determining whether incremental caldera subsidence is dominated by larger, rarer events, or by numerous, repeated, moderate events.
The cumulative subsided volume due to these discrete, meter-scale events, which are encouraged due to Ambrym's small roof aspect ratio (<0.4), is nonnegligible for a caldera of this diameter. The threshold on pressure change or critical volume fraction activating caldera ring-faults is smaller at Ambrym than at calderas with larger roof aspect ratios (Geshi et al., 2014;Geyer et al., 2006). For a caldera with a given diameter, a lower roof aspect ratio (shallower reservoir) translates to a ring-fault with a smaller surface area. This lowers the reservoir depressurization needed for caldera ring-fault reactivation. We note that the estimated depressurization and critical volume fraction depend on the physical properties of the fault, such as the friction coefficient, which are not discussed in detail in this study (Holohan et al., 2015).
It is possible that caldera development by ring-faulting during moderate-sized eruptions occurs at other broad and shallow calderas, such as at Sierra Negra in the Galápagos (Munro & Rowland, 1996). Previous studies have concluded that Sierra Negra's reservoir may be a thin ellipsoid 0.4 km tall and 7 km wide at ∼2-3 km depth, resulting in a roof aspect ratio of ∼0.25-0.4, similar to Ambrym (Amelung et al., 2000;Bell et al., 2021;S. Yun et al., 2006). Much attention has been given to the uplift episodes that cumulate in trap-door faulting along the intracaldera sinuous ridge-fault system in the southwestern portion of the caldera. Little is known, however, about the caldera formation itself, which has been classified as a coherent "piston-type" collapse, but is undated (Jónsson, 2009). Sierra Negra has hosted meter-scale coeruptive subsidence events in both 2005 and 2018 (Abe et al., 2019;Bell et al., 2021;Geist et al., 2008). The magnitude of recent coeruptive subsidence events-5 and 8 m, respectively-makes them difficult to study with InSAR (Casu et al., 2011;S.-H. Yun et al., 2007). During coeruptive subsidence events, fault slip has occurred along portions of the sinuous ridge, but there has been no evidence of reactivation of the bounding caldera ringfaults (Bell et al., 2021).
Future studies concerning subsidence and ring-fault reactivation during moderate-sized eruptions at broad mafic calderas, for example, in the Galápagos, Hawaii, or Iceland are needed to understand the role played by the depth, geometry, and size of the underlying magmatic plumbing system. An outstanding question is whether or not cumulative displacement over geologically significant time scales is net negative or positive. At Ambrym, over the short period of study (20 years), subsidence dominates caldera floor motion (cumulative maximum subsidence of >4 m, including subsidence in 2018-2019 and posteruptive subsidence in 2015, Figures S14 and S15). In addition, geodetic measurements have observed only normal caldera ring-faulting associated with subsidence. Volcanoes, such as those in the Galápagos, can have cycles of preeruptive inflation and coeruptive deflation, which may result in net uplift (i.e., caldera resurgence) or no net displacement over longer periods of time (Bell et al., 2021;Galetto et al., 2019;La Femina et al., 2019).
Small episodic collapses may contribute to long-term geomorphological caldera development if the amount of inelastic subsidence produced by caldera ring-fault reactivation is greater than any inelastic uplift caused by normal faulting produced during intercollapse periods (due to magma replenishment, volatile exsolution, etc.). Ground subsidence at Ambrym may dominate, over decadal time scales, as shown in the uplift/ subsidence ratio measured in the years prior to the December 2018 rift zone intrusion (Shreve, 2020;Shreve et al., 2019). This rift zone intrusion drained a central storage region, causing ∼3 m of subsidence during the event, yet no uplift was measured in the 5 years preceding this event (see cumulative reservoir volume change modeled from geodetic observations spanning 2004-2020, Figure S15). The lack of preeruptive uplift may be due to the open system (i.e., lava lakes persistence and passive degassing), which prevents significant pressurization of the reservoir. Magma replenishment at depth may have occurred in the years before the rift intrusion, as noted by an increase in thermal output (Coppola et al., 2013). However, because the system was open, replenishment may have produced the observed formation of new lava lakes, or the increase in convection vigor and degassing, as opposed to significant reservoir pressurization producing measurable ground displacement (Shreve, 2020). It is possible that Ambrym is an exceptional case where cumulative subsidence results in caldera deepening over hundreds of events due to the lack of intercollapse uplift; however, more extensive monitoring will be needed to confirm this hypothesis. In addition, Ambrym's caldera formed over thousands of years, a time scale that cannot be explored with geodetic observations. Further investigation of this hypothesis is recommended at other calderas that may have formed by similar mechanisms but have a more detailed record of historical eruptions or ring-fault reactivation.

Conclusion
The February 2015 dike intrusion, which resulted in more than 1 m of surface uplift, was concomitant with the partial reactivation of Ambrym's caldera ring-faults, inducing localized subsidence along the caldera rim. Using the 3D Mixed BEM, we conclude that the stress transfer from a depressurizing reservoir most likely promoted ring-faulting (see Figure 6). Incorporating a depressurizing reservoir both decreases the AIC, as well as decreases the external shear stress imposed on the caldera ring-fault in the inversion. The presence of a reservoir feeding the dike intrusion contrasts with previous studies, which neglect this source in their analytical modeling (Hamling & Kilgour, 2020). Our study emphasizes the importance of considering mechanical interactions when sources are in proximity, allowing us to take into account the stress perturbation causing fault slip. This is consistent with studies that systematically investigate errors produced when inverting deformation data from adjacent interacting sources using analytical models (Pascal et al., 2013). Although Pascal et al. (2013) only discusses the interaction between magma reservoirs and dikes, we have shown that source interactions should also be considered for other geological features, such as faults. In addition to considering source interactions, we also compare the volume change in the dike and the reservoir, estimating a magma compressibility at Ambrym that ranges from 4.5 × 10 −10 to 22 × 10 −10 Pa −1 . The critical fraction of extracted material to reactivate ring-faults at Ambrym has an upper SHREVE ET AL.
10.1029/2020JB020277 21 of 26 Figure 6. A simplified schematic illustration summarizing the findings from this study. In 2015, a reservoir (blue spheroid) at depths of 4-5 km b.s.l. fed a dike intrusion and intracaldera fissure eruption. An extraction of less than 7% of magma from the reservoir is sufficient to unclamp the western portion of Ambrym's caldera ring-fault (light green rectangle). During the 2015 eruption, there was an average of 44 cm of slip on the fault plane. We hypothesize that the fault was prestressed by a previous event associated with a depressurizing reservoir (due to an eruption, degassing, etc.), and the reservoir deflation in 2015 brought the fault to failure. Using geodetic models and independent constraints, we conclude that the reservoir feeding the eruption has a minimum volume of 1 km 3 . We cannot put an upper bound of the volume of magma stored in this system (possible lateral extension of magma-storing region outlined with gray dotted lines). Vertical exaggeration of ∼1.5x.
bound of f = 7% for a reservoir with a minimum volume of 1 km 3 . However, we acknowledge that a number of parameters influence f, including uncertainties on 1. the exact timing of onset of ring-fault slip, 2. the amount of stress change on the ring-fault from prior eruptions, reservoir processes, and/or earthquakes, 3. the frictional parameters of the fault, 4. the geometry of the system, especially the connection between the fault and the reservoir, and 5. the reservoir volume's unknown upper bound.
Continued geodetic monitoring of Ambrym, in conjunction with widespread observations at basaltic calderas around the world, will allow us to further constrain how moderate-sized eruptions contribute to basaltic caldera development.

Acknowledgments
We thank the Italian Space Agency (ASI Science project no. 702) and the Japan Aerospace Exploration Agency (JAXA 6th Research agreement no. 3245) for providing access to the radar imagery used in this study. A portion of the ALOS-2 data was provided under a cooperative research contract between Geospatial Information Authority of Japan and JAXA. The ownership of ALOS-2 data belongs to JAXA. The Pléiades data were provided through Dinamis Project 2018-147-Sci, ©CNES (2014©CNES ( , 2015©CNES ( , 2017 and Airbus DS (2017), all rights reserved. We thank Francisco Delgado and Yves Moussallam for discussions which contributed to developing the ideas within this manuscript. Several calculations used the S-CAPAD cluster of IPGP. D. Smittarello, V. Cayol, and V. Pinel acknowledge the support of project ANR-18-CE92-0037. M. Boichu kindly thanks Lieven Clarisse for providing the IASI SO 2 height data and acknowledges the support of the VOLCPLUME ANR project (ANR-15-CE04-0003-01). This project has also received funding from the European Union's Horizon 2020 research and innovation program under the Marie Skłodowska-Curie grant agreement no. 665850. The authors also thank Yosuke Aoki, one anonymous reviewer, and the Associate Editor for their constructive evaluation of this paper. This is IPGP contribution number 4216.