Tara Shreve
email: tshreve@carnegiescience.edu

R Grandin

D Smittarello

V Cayol

V Pinel

M Boichu

Y Morishita

What triggers caldera ring-fault subsidence at Ambrym volcano? Insights from the 2015 dike intrusion and eruption

Ground displacement at Ambrym in February 2015 was caused by a dike intrusion, deflating reservoir, and normal slip on a caldera ring-fault.

• Extracting at most 7% of the magma from Ambrym's reservoir suffices to reactivate the caldera ring-faults.

• Normal slip along Ambrym's ring-fault can occur during moderate-sized eruptions, resulting in subsidence and further caldera development.

Introduction

Although all calderas host collapse structures, instrumentally recorded caldera ringfault activation is rare. A handful of extreme cases of caldera ring-fault activation-catastrophic caldera collapse-have been observed historically [START_REF] Neal | The 2018 rift eruption and summit collapse of Kīlauea volcano[END_REF][START_REF] Gudmundsson | Gradual caldera collapse at Bárdarbunga volcano, Iceland, regulated by lateral magma outflow[END_REF][START_REF] Peltier | Formation of the April 2007 caldera collapse at Piton de La Fournaise volcano: Insights from GPS data[END_REF][START_REF] Geshi | Caldera collapse during the 2000 eruption of Miyakejima Volcano, Japan[END_REF]. Many of these spectacular events included noncoherent, or piecemeal, collapse, where upwards propagating faults break up the collapsing piston (e.g., Bárðarbunga (Ágústsdóttir et al., 2019), Kīlauea [START_REF] Neal | The 2018 rift eruption and summit collapse of Kīlauea volcano[END_REF] K. R. [START_REF] Anderson | Magma reservoir failure and the onset of caldera collapse at Kīlkǒlcano in 2018[END_REF][START_REF] Segall | Caldera collapse geometry revealed by near-field GPS displacements at Kīlauea volcano in 2018[END_REF], Piton de la Fournaise [START_REF] Peltier | Formation of the April 2007 caldera collapse at Piton de La Fournaise volcano: Insights from GPS data[END_REF]).

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 [START_REF] Roche | Experimental study of caldera formation[END_REF]. [START_REF] Walker | Downsag calderas, ring faults, caldera sizes, and incremental caldera growth[END_REF] 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 [START_REF] Druitt | Lithic breccia and ignimbrite erupted during the collapse of Crater Lake Caldera, Oregon[END_REF])) in the literature, despite voluminous historical eruptions [START_REF] Reynolds | Physical volcanology and structural development of Sierra Negra volcano, Isabela Island, Galápagos archipelago[END_REF][START_REF] Munro | Caldera morphology in the western Galápagos and implications for volcano eruptive behavior and mechanisms of caldera formation[END_REF].

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 timescales 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 improve the chances of measuring caldera subsidence and ring-fault activation during moderate-sized eruptions (VEI<3) [START_REF] Pinel | Volcanology: Lessons learned from Synthetic Aperture Radar imagery[END_REF][START_REF] Shreve | From prodigious volcanic degassing to caldera subsidence and quiescence at Ambrym (Vanuatu): the influence of regional tectonics[END_REF].

Caldera ring-faults can be activated when the reservoir pressure drops beneath lithostatic (i.e., is underpressurized) during magma extraction [START_REF] Druitt | On the formation of calderas during ignimbrite eruptions[END_REF]. At basaltic caldera-rift systems, this extraction occurs primarily through lateral dike intrusions [START_REF] Sigmundsson | Calderas collapse as magma flows into rifts[END_REF]. 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 [START_REF] Grandin | September 2005 Manda hararo-dabbahu rifting event, Afar (Ethiopia): Constraints provided by geodetic data[END_REF]. In addition, other factors such as host rock or magma compressibility can reduce the amplitude of the geodetic signal [START_REF] Rivalta | Magma compressibility and the missing source for some dike intrusions[END_REF]. 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 meter 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 Boundary Element Method 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.

Geological Setting

Ambrym Volcano, Vanuatu

Ambrym is a basaltic volcanic island located in the central portion of the New Hebrides subduction zone in Vanuatu (Fig. 1a). The island hosts two rift zones oriented N105 • S, which radiate bilaterally from the island's central 12-km wide caldera [START_REF] Mccall | The Geology and Geophysics of the Ambrym Caldera, New Hebrides[END_REF] (Fig. 1b). The caldera has been previously interpreted as resulting from the collapse of a giant tuff cone during a sequence of explosive phreatomagmatic eruptions [START_REF] Robin | Giant tuff cone and 12-km-wide associated caldera at Ambrym Volcano (Vanuatu, New Hebrides Arc)[END_REF]. This view was later challenged by [START_REF] Cronin | Where are the giant tuff cone and ignimbrites of Ambrym? A more conventional story of mafic volcanism at Ambrym Volcano, Vanutu[END_REF], 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 travelling more than 20 km into the SE rift zone, resulting in a submarine eruption that emitted basaltic pumice onto the nearby shoreline [START_REF] Shreve | From prodigious volcanic degassing to caldera subsidence and quiescence at Ambrym (Vanuatu): the influence of regional tectonics[END_REF]. Coeval with the 2018 rift zone intrusion, the caldera floor subsided by more than 2 meters. The caldera ring-faults also reactivated, most notably in the northern half of the caldera. A depressurizing deformation source was modelled at between 1-5 km beneath the summit, depending on how many sources were included [START_REF] Hamling | Large-scale drainage of a complex magmatic system: Observations from the 2018 eruption of Ambrym volcano, Vanuatu[END_REF][START_REF] Shreve | From prodigious volcanic degassing to caldera subsidence and quiescence at Ambrym (Vanuatu): the influence of regional tectonics[END_REF]. 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 2000 years ago [START_REF] Robin | Giant tuff cone and 12-km-wide associated caldera at Ambrym Volcano (Vanuatu, New Hebrides Arc)[END_REF]. Due to these observations of meter-scale caldera subsidence, in conjunction with historical documentation of reoccurring rift zone intrusions, recent studies [START_REF] Hamling | Large-scale drainage of a complex magmatic system: Observations from the 2018 eruption of Ambrym volcano, Vanuatu[END_REF][START_REF] Shreve | From prodigious volcanic degassing to caldera subsidence and quiescence at Ambrym (Vanuatu): the influence of regional tectonics[END_REF] 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 moderatesized eruptions, such as an intra-caldera 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) [START_REF] Németh | Volcanic craters, pit craters and high-level magma-feeding systems of a mafic island-arc volcano: Ambrym, Vanuatu, South Pacific[END_REF]. In addition to lava lake activity, occasional intra-caldera fissure eruptions have occurred, most notably in 1986 and 1988˘-89 [START_REF] Eissen | Chronology of the historic volcanic activity of the New Hebrides island arc from 1595 to 1991[END_REF] (Fig. 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 [START_REF] Robin | Giant tuff cone and 12-km-wide associated caldera at Ambrym Volcano (Vanuatu, New Hebrides Arc)[END_REF][START_REF] Eissen | Chronology of the historic volcanic activity of the New Hebrides island arc from 1595 to 1991[END_REF].

February 2015 Eruption

On 20-21 February 2015, another intra-caldera fissure eruption occurred, with a main fissure located ∼3 km SE of Marum (Fig. 1d) (Coppola, Laiolo, & Cigolini, 2016;[START_REF] Hamling | Goldilocks conditions required for earthquakes to trigger basaltic eruptions: Evidence from the 2015 Ambrym eruption[END_REF]. 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 19 February 2015 at 13:18 UTC with a hypocenter 30 km southeast of Ambyrm's craters and a depth of ∼5 km (Fig. 1a, 16.50 • S, 168.28 • E, according to the location from the Oceania Regional Seismic NETwork (ORSNET), event ID ird2015dmnz).

The earthquake had a reverse focal mechanism and moment tensor solution that included only 61% double couple component (Fig. 1a). Thermal anomalies were detected by both MODVOLC (R. Wright et al., 2004;[START_REF] Wright | MODVOLC: 14 years of autonomous observations of effusive volcanism from space[END_REF] and the Middle Infrared Observation of Volcanic Activity (MIROVA) system [START_REF] Coppola | Enhanced volcanic hot-spot detection using MODIS IR data: Results from the MIROVA system[END_REF]. According to Coppola, Laiolo, and Cigolini (2016), 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 (Coppola, Laiolo, & Cigolini, 2016) (see Supporting Information Text S1 for confirmation of this onset time using the Hybrid Single-Particle Lagrangian Integrated Trajectory model (HYSPLIT) [START_REF] Stein | NOAA's hysplit atmospheric transport and dispersion modeling system[END_REF][START_REF] Crawford | Initializing HYSPLIT with satellite observations of volcanic ash: A case study of the 2008 Kasatochi eruption[END_REF] and information on the SO 2 plume height derived from Infrared Atmospheric Sounding Interferometer (IASI) observations [START_REF] Clarisse | The 2011 Nabro eruption, a SO 2 plume height analysis using IASI measurements[END_REF]). Coppola, Laiolo, and Cigolini (2016) 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, travelling ∼3 km.

To directly constrain the lava flow volume, we calculate a digital elevation model (DEM) using MicMac software [START_REF] Rupnik | 3D reconstruction from multi-view VHR-satellite images in MicMac[END_REF] and post-eruption 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 (Fig. 1d, Supporting Information Fig. S1, Text S2 and Table S1). Uncertainties were calculated using the methods of [START_REF] Bagnardi | High-resolution digital elevation model from tri-stereo Pleiades-1 satellite imagery for lava flow volume estimates at Fogo Volcano[END_REF] and [START_REF] Favalli | Evolution of an active lava flow field using a multitemporal LIDAR acquisition[END_REF]. This value is ∼2 times larger than the value estimated by Coppola, Laiolo, and Cigolini (2016).

We use the Ozone Mapping and Profiler Suite (OMPS) observations to conservatively estimate ∼40 kt of sulfur dioxide (SO 2 ) emitted during the eruption (Fig. 1e, Supporting Information Text S3). For this SO 2 mass estimation, we used the OMPS TRM product, which assumes a center of mass altitude (CMA) of ∼7.5 km [START_REF] Li | New-generation NASA Aura Ozone Monitoring Instrument (OMI) volcanic SO 2 dataset: Algorithm description, initial results, and continuation with the Suomi-NPP Ozone Mapping and Profiler Suite (OMPS)[END_REF]. 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 25 February, SO 2 degassing returned to passive background levels (∼7 kt/day, according to [START_REF] Carn | A decade of global volcanic SO 2 emissions measured from space[END_REF]). 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 synthetic aperture radar (SAR) datasets to measure the co-eruptive ground displacement (Table S2). A joint analysis of differential In-SAR and Multiple Aperture Interferometry from ALOS-2, and pixel offset tracking from COSMO-SkyMed (CSK), allows us to decompose the 3D co-eruptive displacement field.

We use the 3D displacements to interpret the deformation signal, determining the possible location and nature of the active deformation sources.

Processing Methods

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 24 January to 21 March 2015 with the Interferometric SAR scientific computing environment (ISCE) [START_REF] Rosen | The InSAR scientific computing environment[END_REF]. Filtering and unwrapping are performed with NSBAS modules [START_REF] Doin | Presentation of the small baseline NSBAS processing chain on a case example: the Etna deformation monitoring from 2003 to 2010 using Envisat data[END_REF][START_REF] Rosen | Updated Repeat Orbit Interferometry Package Released[END_REF][START_REF] Grandin | Long-term growth of the Himalaya inferred from interseismic InSAR measurement[END_REF]. After geocoding, the final pixel posting is ∼14 m×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 8 February 2015, due to inadequate (less than 50%) burst synchronization on the ground for acquisitions before this date [START_REF] Lindsey | Line-of-sight displacement from ALOS-2 interferometry: M w 7.8 Gorkha Earthquake and M w 7.3 aftershock[END_REF][START_REF] Natsuaki | SAR interferometry using ALOS-2 PALSAR-2 data for the M w 7.8 Gorkha, Nepal earthquake[END_REF]. This was due to an issue with ALOS-2's navigation system, which was fixed on 8 February 2015 [START_REF] Lindsey | Line-of-sight displacement from ALOS-2 interferometry: M w 7.8 Gorkha Earthquake and M w 7.3 aftershock[END_REF]. As a result, standard ScanSAR-to-ScanSAR processing is replaced by ScanSAR-to-stripmap processing for descending images spanning 17 January 2015 to 14 March 2015 [START_REF] Ortiz | ScanSAR-to-Stripmap Mode Interferometry Processing Using ENVISAT / ASAR Data[END_REF][START_REF] Natsuaki | SAR interferometry using ALOS-2 PALSAR-2 data for the M w 7.8 Gorkha, Nepal earthquake[END_REF][START_REF] Kobayashi | Detailed crustal deformation and fault rupture of the 2015 Gorkha earthquake, Nepal, revealed from ScanSARbased interferograms of ALOS-2[END_REF]. This is implemented with the GSISAR software developed by the Geospatial Information Authority of Japan [START_REF] Kobayashi | Detailed crustal deformation and fault rupture of the 2015 Gorkha earthquake, Nepal, revealed from ScanSARbased interferograms of ALOS-2[END_REF][START_REF] Tobita | Development of SAR interferometry analysis and its application to crustal deformation study[END_REF]. A 30-m mesh Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) Global Digital Elevation Model (GDEM) was used to remove the topographic fringes [START_REF] Tachikawa | Characteristics of ASTER GDEM Version 2[END_REF]. The post-interferogram formation steps (multilooking, filtering, unwrapping and geocoding) are performed with NS-BAS modules, as described in the Supporting Information Text S4. The descending interferogram is multilooked 1 time in range and 8 times in azimuth, resulting in a final pixel posting of ∼30 m×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) [START_REF] Bechor | Measuring two-dimensional movements using a single InSAR pair[END_REF][START_REF] Liang | Measuring Azimuth Deformation with L-Band ALOS-2 ScanSAR Interferometry[END_REF], in order to derive the along-track, or azimuth, displacement using ascending SM3 images spanning 24 January to 4 April 2015. These measurements are then used in the 3D displacement decomposition. Azimuth common-band filtering with a normalized squint of 0.66 generates sub-aperture single-look-complex (SLC) 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 post-processing steps include filtering [START_REF] Goldstein | Radar Interferogram Filtering for Geophysical Applications[END_REF], unwrapping, and geocoding with a pixel spacing of ∼85 m×70 m.

Pixel Offset Tracking

Sub-pixel offset tracking measures surface displacement by finding the cross-correlation peak of two SAR image amplitudes [START_REF] Michel | Measuring ground displacements from SAR amplitude images : application to the Landers earthquake[END_REF]. Pixel offsets are less accurate than InSAR [START_REF] Michel | Measuring ground displacements from SAR amplitude images : application to the Landers earthquake[END_REF]. 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 to 25 February 2015, we calculate pixel offsets in both range and azimuth. We also calculate pixel offsets from the ALOS-2 ascending pair spanning 24 January to 21 March. 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 post-processing 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. Wright et al. (2004), we invert for the 3D co-eruptive 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 dataset was weighted equally in the inversion, which assumes similar accuracy for all measurements. InSAR has a higher accuracy than the other datasets, but the equal weighting ensures that contributions from the pixel offsets and MAI, such as displacement along the satellite azimuth, are not under-prioritized in the 3D decomposition [START_REF] Grandin | Significance of near-field geodetic data and diversified geometries of acquisition to model coseismic slip and refine shallow slip deficit: the case of the Balochistan earthquake[END_REF]. Notable characteristics of the 3D displacement field allow us to hypothesize which deformation sources were active during this event.

Co-eruptive Displacement Field Description

The ascending and descending interferograms measure, respectively, up to 1.5 and 1 meters of line-of-sight (LOS) shortening (ground motion towards 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 co-eruptive 3D displacement field:

1. There exists a N-S asymmetry and discontinuity across the eruptive fissures in all components (see Figure 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 towards 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 21 February, during the onset of the eruption (Fig. 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 (2 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 modelled by a shallow (4-5 km depth) depressurizing sill [START_REF] Shreve | From prodigious volcanic degassing to caldera subsidence and quiescence at Ambrym (Vanuatu): the influence of regional tectonics[END_REF].

Inversion Scheme

The displacement field's complexity hints at contributions from multiple sources.

Geodetic modelling of volcanic deformation often capitalizes on computationally efficient analytical solutions to describe deformation due to pressure sources [START_REF] Mogi | Relations between the eruptions of various volcanoes and the deformations of the ground surfaces around them[END_REF], as well as shear or tensile dislocations [START_REF] Okada | Surface Deformation Due to Shear and Tensile Faults in a Half-Space[END_REF]. 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 (pre-existing weaknesses and faults due to caldera collapse [START_REF] Acocella | Understanding caldera structure and development: An overview of analogue models compared to natural calderas[END_REF])), we choose to proceed with the 3D Mixed Boundary Element Method (BEM) (see Section 4.2 for more details) [START_REF] Cayol | 3D mixed boundary elements for elastostatic deformation field analysis[END_REF][START_REF] Sambridge | Exploring multidimensional landscapes without a map[END_REF]. 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 modelled sources. By combining the Mixed BEM forward models with a non-linear 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 modelled 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 signal-to-noise ratio. However, we project modelled surface displacements into the LOS for all available datasets a posteriori (Fig. S4). The data is subsampled using an adaptive quadtree decomposition algorithm [START_REF] Walstead | Fractal and wavelet image compression techniques[END_REF][START_REF] Jónsson | Fault slip distribution of the 1999 M w 7.1 Hector Mine, California, earthquake, estimated from satellite radar and GPS measurements[END_REF], as described in the Supporting Information Text S5, Table S4 andFig. S5.

Displacement values at subsampled points may be spatially correlated due to phase contributions from atmospheric, ionospheric, or other noise sources [START_REF] Sudhaus | Improved source modelling through combined use of InSAR and GPS under consideration of correlated data errors: Application to the June 2000 Kleifarvatn earthquake, Iceland[END_REF]. 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

C d (k 1 , k 2 ) = σ 2 d • exp(-||k1,k2|| ac
), 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 ) [START_REF] Tarantola | Inverse problem theory: Method for data fitting and model parameter estimation[END_REF]. Following [START_REF] Fukushima | Finding realistic dike models from interferometric synthetic aperture radar data: The February 2000 eruption at Piton de la Fournaise[END_REF]; [START_REF] Smittarello | Magma Propagation at Piton de la Fournaise From Joint Inversion of InSAR and GNSS[END_REF], we estimated σ 2 d = 3 × 10 -4 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 [START_REF] Cayol | 3D mixed boundary elements for elastostatic deformation field analysis[END_REF][START_REF] Sambridge | Exploring multidimensional landscapes without a map[END_REF].

In the 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 TanDEM-X DEM (Fig. S6a) (B. [START_REF] Wessel | TanDEM-X Ground Segment -DEM Products Specification Document[END_REF]. 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 Mixed BEM combines the displacement discontinuity and direct displacement methods [START_REF] Lachat | Effective numerical treatment of boundary integral equations: A formulation for three-dimensional elastostatics[END_REF][START_REF] Sokolnikoff | Mathematical theory of elasticity[END_REF][START_REF] Crouch | Solution of plane elasticity problems by the displacement discontinuity method[END_REF]. A description of these methods can be found in the Supporting Information Text S6. Nine parameters control the fracture geometry (Fig. S7) [START_REF] Fukushima | Evolution of magma conduits during the 1998-2000 eruptions of Piton de la Fournaise volcano, Réunion Island[END_REF]. 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 ringfault, 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 tradeoff between accuracy and computational efficiency of the forward model, controlled by the number of meshed elements used to represent the boundaries (see Fig. S6b for a plot of model run duration versus 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 (Fig. 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 [START_REF] Cayol | Three-dimensional modeling of the 1983-1984 eruption at Piton de la Fournaise Volcano, Réunion Island[END_REF] and laboratory measurements [START_REF] Heap | Towards more realistic values of elastic moduli for volcano modelling[END_REF]. The influence of the chosen shear modulus on the estimated pressure change and reservoir volume will be discussed further in Section 6.2.3.

Non-linear 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 [START_REF] Sambridge | Geophysical inversion with a neighbourhood algorithm-I. Searching a parameter space[END_REF][START_REF] Fukushima | Finding realistic dike models from interferometric synthetic aperture radar data: The February 2000 eruption at Piton de la Fournaise[END_REF]. The neighborhood algorithm utilises 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 [START_REF] Sambridge | Exploring multidimensional landscapes without a map[END_REF][START_REF] Tridon | Inversion of coeval shear and normal stress of Piton de la Fournaise flank displacement[END_REF]. Each subsequent iteration picks N 2 combinations of parameters, and the N 2 computations can be performed in parallel [START_REF] Fukushima | Finding realistic dike models from interferometric synthetic aperture radar data: The February 2000 eruption at Piton de la Fournaise[END_REF][START_REF] Smittarello | Magma Propagation at Piton de la Fournaise From Joint Inversion of InSAR and GNSS[END_REF]. We take N 2 = 48 to optimize this parallelization, which was performed on a cluster using 3 nodes, with 16 CPUs per node. A particular iteration searches within the N 2 multi-dimensional Voronoi cells with the lowest misfit from the previous iteration.

The misfit is defined as

χ 2 = (u o -u m ) T C -1 d (u o -u m )
, where u o is the data, u m is the modelled displacements, and C d is the data covariance matrix, as defined in Section 4.1. This iterative procedure continues until a threshold criterion, σ N LST , is reached, based on the standard deviation of the misfit values from the last N LST forward calculations. We set N LST = N 2 . To calculate posterior probability density functions (PPD's), marginal PPD's, mean model parameters, and model uncertainties, we follow the bayesian inference framework [START_REF] Tarantola | Inverse problem theory: Method for data fitting and model parameter estimation[END_REF] described in Sambridge (1999a) and implemented by [START_REF] Fukushima | Finding realistic dike models from interferometric synthetic aperture radar data: The February 2000 eruption at Piton de la Fournaise[END_REF] (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 [START_REF] Zebker | On the derivation of coseismic displacement fields using differential radar interferometry: the Landers earthquake[END_REF]. Therefore, we mask the displacement signal, fit a ramp to the non-deforming areas, and subtract the ramp from the ascending interferogram using General Mapping Tools (P. [START_REF] Wessel | The Generic Mapping Tools version 6[END_REF]. 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 non-linearly 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 [START_REF] Akaike | A new look at the statistical model identification[END_REF]. In each inversion, we use the same subsampled datasets and data covariance matrix from Section 4.1. Therefore, to compare the AIC of inversions, we can simply compare AIC = 2k + χ 2 .

Inversion Results

Dike (Inversion 1)

According to (1) in Section 3.3, the N-S asymmetry across the eruptive fissure implies a dike dipping towards the SW. The first inversion estimates 10 parameters nonlinearly-the 9 dike geometry parameters and the pressure change. Table 2 lists all inverted, fixed, and post-processed parameters. The best-fit model, as well as data, synthetics and residuals, are shown in Figure 3. After the inversion is complete, the bestfit 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 [START_REF] Fukushima | Finding realistic dike models from interferometric synthetic aperture radar data: The February 2000 eruption at Piton de la Fournaise[END_REF]. 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 PPD's, are listed in the Supporting Information for all inversions (Table S5 andFigures S9, S10, andS11). 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 (Fig. 4a, 3a, Table S7).

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) in Section 3.3. The second inversion includes both a dike, as well as a frictionless fault, whose top edge is predefined and intersects the topography mesh. We invert for 3 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 non-linearly 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 Fig. S7 for a description of these parameters). This simplification is chosen because analogue and numerical models have found caldera ring-faults to be vertical at depth [START_REF] Beauducel | [END_REF][START_REF] Gudmundsson | Gradual caldera collapse at Bárdarbunga volcano, Iceland, regulated by lateral magma outflow[END_REF][START_REF] Liu | Influence of ring faulting in localizing surface deformation at subsiding calderas[END_REF]. We also only include the portion of the fault that we hypothesize was reactivated in 2015 (See Section 6.1). Because the fault and dike orientations are perpendicular, pressurization of the dike results in closing (fracture wall interpen-etration) on some sections of the fault plane (Fig. S8). We impose a nonnegativity constraint to avoid fracture wall interpenetration. This nonnegativity constraint is enforced in the Mixed BEM through Lagrange multipliers, as described in [START_REF] Cayol | Sheared sheet intrusions as mechanism for lateral flank displacement on basaltic volcanoes: Applications to Réunion Island volcanoes[END_REF] (from hereon we call this constraint "preventing interpenetration"). As a result, computation time of a single forward model increases from 30 seconds 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 [START_REF]-D modelling of Campi Flegrei ground deformations: Role of Caldera boundary discontinuities[END_REF] in Section 3.3).

Therefore, a final inversion includes a thin (semi-minor 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 post-intrusion deflation source estimated in [START_REF] Shreve | From prodigious volcanic degassing to caldera subsidence and quiescence at Ambrym (Vanuatu): the influence of regional tectonics[END_REF]. Inverting for these parameters would result in an unreasonably long inverison. As shown by [START_REF] Sambridge | Exploring multidimensional landscapes without a map[END_REF], 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 datasets (e.g., [START_REF] Shreve | From prodigious volcanic degassing to caldera subsidence and quiescence at Ambrym (Vanuatu): the influence of regional tectonics[END_REF]; [START_REF] Hamling | Large-scale drainage of a complex magmatic system: Observations from the 2018 eruption of Ambrym volcano, Vanuatu[END_REF]; [START_REF] Hamling | Goldilocks conditions required for earthquakes to trigger basaltic eruptions: Evidence from the 2015 Ambrym eruption[END_REF]; [START_REF] Allard | Prodigious emission rates and magma degassing budget ofmajor, trace and radioactive volatile species from Ambrym basaltic volcano, Vanuatu island Arc[END_REF]; [START_REF] Legrand | Interpretation of very long period tremors at Ambrym volcano, Vanuatu, as quasi-static displacement field related to two distinct magmatic sources[END_REF] for reservoir depth) and post-processing (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 [START_REF] Amoruso | Shape and volume change of pressurized ellipsoidal cavities from deformation and seismic data[END_REF], resulting in the high ∆P res obtained by the inversion given the thin reservoir. After running Mixed BEM forward models for spheroids with various thicknesses, we find that the difference in both the surface displacement and the source volume change are within reasonable error bounds for all aspect ratios (see Section 6.2.2 and Fig. 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 PPD's 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 = Table 1. The number of free parameters, number of total forward models, AIC (calculated as AIC = 2k + χ 2 , where k is the number of inverted parameters and χ 2 is the misfit). Each inversion ended when the standard deviation criterion (<0.5) was reached. -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 one meter 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 (Fig. 4), due at least in part to the lack of sensitivity of the surface displacement to pressure changes from deep sources [START_REF] Du | Comparison of various inversion techniques as applied to the determination of a geophysical deformation model for the 1983 Borah Peak earthquake[END_REF]. The mean fault dip estimated in both Inversion 2 and 3 are within a few degrees, however, the dip tends towards a shallower value. We decide to constrain the minimum dip to 65 • based on numerical and analogue modelling, 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) (Roche et al. (2001) and references within). According to geodetic modelling from the 2018 eruption [START_REF] Shreve | From prodigious volcanic degassing to caldera subsidence and quiescence at Ambrym (Vanuatu): the influence of regional tectonics[END_REF][START_REF] Hamling | Large-scale drainage of a complex magmatic system: Observations from the 2018 eruption of Ambrym volcano, Vanuatu[END_REF], as well as depths of Ambrym's shallow reservoir estimated with melt inclusion entrapment pressures and very long period tremors [START_REF] Allard | Prodigious emission rates and magma degassing budget ofmajor, trace and radioactive volatile species from Ambrym basaltic volcano, Vanuatu island Arc[END_REF][START_REF] Legrand | Interpretation of very long period tremors at Ambrym volcano, Vanuatu, as quasi-static displacement field related to two distinct magmatic sources[END_REF]), Ambrym's reservoir roof most likely has an R a < 0.4. In extreme cases, caldera sagging can result in material slumping along shallowly-dipping normal detachments formed within the footwall of previously formed normal faults [START_REF] Holohan | Mechanical and geometric controls on the structural evolution of pit crater and caldera subsidence[END_REF], 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 the Table S5 and Table S6.

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 sourcesa 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. In the context of a caldera, to test if stress perturbations due to dike opening promote normal caldera ring-faulting, we calculate the Coulomb stress change on the fault plane [START_REF] King | Static Stress Change and the Triggering of Earthquakes[END_REF][START_REF] Lin | Stress triggering in thrust and subduction earthquakes and stress interaction between the southern San Andreas and nearby thrust and strike-slip faults[END_REF]. The Coulomb stress change is 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 [START_REF] Nostro | Twoway coupling between Vesuvius eruptions and southern Apennine earthquakes, Italy, by elastic stress transfer[END_REF]. 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 [START_REF] Rubin | Dike-induced faulting in rift zones of Iceland and Afar[END_REF]. 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 modelling 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 [START_REF] Buck | Tectonic stress and magma chamber size as controls on dike propagation: Constraints from the 1975-1984 Krafla rifting episode[END_REF][START_REF] Pinel | The effect of edifice load on magma ascent beneath a volcano[END_REF][START_REF] Corbi | How caldera collapse shapes the shallow emplacement and transfer of magma in active volcanoes[END_REF]. 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 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 (Fig. 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 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 MPa 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 non-zero 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 pre-existing 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 pre-existing external stress on the fault. For example, Ambrym's persistently high rates of passive de-gassing over the past decade may have caused reservoir depressurization, if there is no magma replenishment from a deeper source [START_REF] Girona | On depressurization of volcanic magma reservoirs by passive degassing[END_REF]. 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.

Caldera Ring-fault Reactivation

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 [START_REF] Sigmundsson | Calderas collapse as magma flows into rifts[END_REF]. The plumbing system's geometry influences the critical fraction of magma needed to trigger collapse (defined as f crit = -∆q f V , where ∆q f is the magma extraction volume before the onset of collapse and V is the reservoir volume). Analytical, analogue, and numerical models conclude that thresholds are lower for shallower, broader calderas (low roof aspect ratios, R a = h 2r , where h is reservoir depth and r is the reservoir radius) [START_REF] Roche | Experimental study of caldera formation[END_REF][START_REF] Geshi | Evaluating volumes for magma chambers and magma withdrawn for caldera collapse[END_REF][START_REF] Holohan | Mechanical and geometric controls on the structural evolution of pit crater and caldera subsidence[END_REF]. 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 analogue modelling (K. R. [START_REF] Anderson | Magma reservoir failure and the onset of caldera collapse at Kīlkǒlcano in 2018[END_REF] and references within). Ambrym has an especially low roof aspect ratio (<0.4) [START_REF] Shreve | From prodigious volcanic degassing to caldera subsidence and quiescence at Ambrym (Vanuatu): the influence of regional tectonics[END_REF][START_REF] Hamling | Large-scale drainage of a complex magmatic system: Observations from the 2018 eruption of Ambrym volcano, Vanuatu[END_REF], 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 pre-existing caldera ring-faults at a broad, shallow basaltic caldera.

Pressure Change and Minimum Reservoir Volume

Using the Mixed BEM and independent obsevations, 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 semi-axes a, b, and c. Because ∆V res ∝ ∆PresV µ , a tradeoff exists between ∆P res and V (the full relationship for various reservoir geometries can be found in [START_REF] Amoruso | Shape and volume change of pressurized ellipsoidal cavities from deformation and seismic data[END_REF]). Here, c represents the "height" or "thickness" of the reservoir, which we fix in the inversion because of this tradeoff. We explore the range of reasonable V by setting the semi-minor axis c of the reservoir estimated in Inversion 3 to values that range between 2.5 m and 1300 m, adjusting ∆P res to obtain a ∆V res ≈ -15 × 10 6 m 3 (as estimated in Inversion 3). The residuals between the forward Mixed BEM calculations for various aspect ratios and the synthetic displacement field from Inversion 3 are shown in Figure S12. Val-ues 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 [START_REF] Shreve | Crustal deformation at Ambrym (Vanuatu) imaged with satellite geodesy: constraints on magma storage, migration, and outgassing[END_REF] and after Ambrym's 2018 rift zone intrusion [START_REF] Shreve | From prodigious volcanic degassing to caldera subsidence and quiescence at Ambrym (Vanuatu): the influence of regional tectonics[END_REF].

We can put a lower bound on the reservoir's negative pressure change, because we do not observe a full-scale 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 Roche and Druitt (2001) (assuming vertical faults) and a maximum host rock shear strength τ max of 100 MPa [START_REF] Touloukian | Physical properties of rocks and minearls. data serises on material properties[END_REF][START_REF] Schön | Fundamentals and principles of petrophysics, handbook of geophysical exploration, seismic exploration[END_REF], 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 [START_REF] Sparks | Causes and consequences of pressurisation in lava dome eruptions[END_REF][START_REF] Anderson | Bayesian inversion of data from effusive volcanic eruptions using physics-based models: Application to Mount St. Helens 2004-2008[END_REF][START_REF] Delgado | Physicochemical models of effusive rhyolitic eruptions constrained with InSAR and DEM data: A case study of the 2011-2012 Cordón Caulle eruption[END_REF]. Given the estimated reservoir dimensions of a ≈ 1340 m from Inversion 3, to obtain a -∆P < 175 MPa, c can range from 5 to 1300 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 [START_REF] Allard | Prodigious emission rates and magma degassing budget ofmajor, trace and radioactive volatile species from Ambrym basaltic volcano, Vanuatu island Arc[END_REF] estimates Ambrym's minimum reservoir volume to be 0.5 km 3 . Allard et al. (2015) obtains 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, [START_REF] Allard | Prodigious emission rates and magma degassing budget ofmajor, trace and radioactive volatile species from Ambrym basaltic volcano, Vanuatu island Arc[END_REF] assumes 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 [START_REF] Carn | A decade of global volcanic SO 2 emissions measured from space[END_REF] (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 [START_REF] Carn | A decade of global volcanic SO 2 emissions measured from space[END_REF] 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 modelling, 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 [START_REF] Shreve | From prodigious volcanic degassing to caldera subsidence and quiescence at Ambrym (Vanuatu): the influence of regional tectonics[END_REF][START_REF] Hamling | Large-scale drainage of a complex magmatic system: Observations from the 2018 eruption of Ambrym volcano, Vanuatu[END_REF]. 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 of reservoir volume change to the volume of intruded material, r v , we can estimate to the first-order the magma compressibility [START_REF] Rivalta | Magma compressibility and the missing source for some dike intrusions[END_REF][START_REF] Rivalta | Evidence that coupling to magma chambers controls the volume history and velocity of laterally propagating intrusions[END_REF]. 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 [START_REF] Johnson | Comment on 'volume of magma accumulation or withdrawal estimated from surface uplift or subsidence, with application to the 1960 collapse of Kīlauea volcano' by p.t. delaney and d.f. mctigue[END_REF]. [START_REF] Rivalta | Magma compressibility and the missing source for some dike intrusions[END_REF]; [START_REF] Rivalta | Evidence that coupling to magma chambers controls the volume history and velocity of laterally propagating intrusions[END_REF] derived the generic formula for r v based on chamber compressibility β c and magma compressibility β m , such that

r v = ∆V s ∆V c = 1 + β m β c , (1)
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, 2011) interpolated these estimations using finite element calculations for intermediate aspect ratios.

β c = 1 V ∂V ∂P ,
The reservoir that drained during the 2015 Ambrym eruption can be modelled as an ellipsoid with a = b and an aspect ratio c a . We calculate chamber compressibility using analytical expressions from [START_REF] Amoruso | Shape and volume change of pressurized ellipsoidal cavities from deformation and seismic data[END_REF], finite element calculations (K. [START_REF] Anderson | Physics-based models of ground deformation and extrusion rate at effusively erupting volcanoes[END_REF], and the 3D Boundary Element Method. 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

β c = 3 µ a c 1 2π - 1 5 . (2)
Although Equation 2 fails for shallow penny-shaped sills (K. [START_REF] Anderson | Physics-based models of ground deformation and extrusion rate at effusively erupting volcanoes[END_REF], 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 upperbound, chamber compressibility can be approximated by a sphere within 25% error when aspect ratios c a range between ∼0.6-10 (K. [START_REF] Anderson | Physics-based models of ground deformation and extrusion rate at effusively erupting volcanoes[END_REF], such that

β c = 3 4µ . (3)
Using Equation 2 and 3, it follows that β c ranges from 3.8 × 10 -10 -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 Boundary Element Method. β c in this case ranges from 3.7×10 -10 -18×10 -10 . The differences between the analytical and the computational values are consistent with K. [START_REF] Anderson | Physics-based models of ground deformation and extrusion rate at effusively erupting volcanoes[END_REF], 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 r v = (∆V dike +∆V erupted ) ∆Vres ≈ 2.2 (implying β m ≈ β c ).
It follows that magma compressibility ranges from 4.5×10 -10 -22×10 -10 Pa -1 . Although spanning almost one order of magnitude, we can conclude that this value is higher than that for degassed basalts [START_REF] Spera | Physical properties of magma[END_REF]. 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 [START_REF] Hamling | Goldilocks conditions required for earthquakes to trigger basaltic eruptions: Evidence from the 2015 Ambrym eruption[END_REF]. 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 -15×10 -10 Pa -1 , also higher than values estimated assuming gas-poor basalt (K. R. [START_REF] Anderson | Magma reservoir failure and the onset of caldera collapse at Kīlkǒlcano in 2018[END_REF][START_REF] Segall | Caldera collapse geometry revealed by near-field GPS displacements at Kīlauea volcano in 2018[END_REF].

K. R. [START_REF] Anderson | Magma reservoir failure and the onset of caldera collapse at Kīlkǒlcano in 2018[END_REF]; [START_REF] Segall | Caldera collapse geometry revealed by near-field GPS displacements at Kīlauea volcano in 2018[END_REF] 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 R = βm βc . Similarly, at Erta Ale, another volcano which hosts a lava lake, no significant cumulative post-intrusion subsidence was measured after the 2017 intrusion, indicating either highly compressible magma or magma fed from a deep source [START_REF] Moore | The 2017 Eruption of Erta 'Ale Volcano, Ethiopia: Insights into the Shallow Axial Plumbing System of an Incipient Mid-Ocean Ridge[END_REF]. 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 [START_REF] Allard | Prodigious emission rates and magma degassing budget ofmajor, trace and radioactive volatile species from Ambrym basaltic volcano, Vanuatu island Arc[END_REF], who concludes that Ambrym's degassing occurs in a "closed system" [START_REF] Edmonds | New geochemical insights into volcanic degassing[END_REF]. "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 [START_REF] Shreve | From prodigious volcanic degassing to caldera subsidence and quiescence at Ambrym (Vanuatu): the influence of regional tectonics[END_REF]. This balance between the degassed sulfur mass and erupted lava volume has already been documented at other effusive volcanic eruptions [START_REF] Barnie | A multidisciplinary study of the final episode of the Manda Hararo dyke sequence, Ethiopia, and implications for trends in volcanism during the rifting cycle[END_REF]), yet this is not necessarily the case in explosive eruptions, as emphasized by the "excess" sulfur problem [START_REF] Wallace | Volcanic SO 2 emissions and the abundance and distribution of exsolved gas in magma bodies[END_REF][START_REF] Kilbride | Observing eruptions of gasrich compressible magmas from space[END_REF].

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 [START_REF] Heap | Towards more realistic values of elastic moduli for volcano modelling[END_REF]. 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 [START_REF] Heap | Towards more realistic values of elastic moduli for volcano modelling[END_REF]. 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 mod-ulus 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 [START_REF] Heap | Towards more realistic values of elastic moduli for volcano modelling[END_REF], or with seismic wave velocities [START_REF] Ellis | A future magma inflation event under the rhyolitic Taupo volcano, New Zealand: Numerical models based on constraints from geochemical, geological, and geophysical data[END_REF][START_REF] Grandin | Sequence of rifting in Afar, Manda-Hararo rift, Ethiopia, 2005-2009: Time-space evolution and interactions between dikes from interferometric synthetic aperture radar and static stress change modeling[END_REF][START_REF] Masterlark | Nonlinear estimation of geometric parameters in FEMs of volcano deformation: Integrating tomography models and geodetic data for Okmok volcano, Alaska[END_REF][START_REF] Wauthier | Magma pathways and their interactions inferred from InSAR and stress modeling at Nyamulagira Volcano[END_REF][START_REF] Albino | The Role of Pore Fluid Pressure on the Failure of Magma Reservoirs: Insights From Indonesian and Aleutian Arc Volcanoes[END_REF]). 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 post-eruptive uplift related to magmatic recharge was measured, despite SAR acquisitions every 6 days [START_REF] Shreve | From prodigious volcanic degassing to caldera subsidence and quiescence at Ambrym (Vanuatu): the influence of regional tectonics[END_REF]. In addition, the caldera floor subsided in the months following the 2015 eruption, as measured with both ALOS-2 and Sentinel-1 (Fig. S15 andS16). 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 ringfault 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 [START_REF] Shreve | From prodigious volcanic degassing to caldera subsidence and quiescence at Ambrym (Vanuatu): the influence of regional tectonics[END_REF][START_REF] Hamling | Large-scale drainage of a complex magmatic system: Observations from the 2018 eruption of Ambrym volcano, Vanuatu[END_REF].

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 1000×10 6 m 3 . Because the reservoir deflation and caldera ring-fault reactivation occurred before the end of the eruption (Fig. 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 V crit = -∆V V is thus at most 3%, and, after taking into account the compressibility of the system, the critical fraction f = -∆q 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 analogue models. [START_REF] Geyer | Relationship between caldera collapse and magma chamber withdrawal : An experimental approach[END_REF] estimates 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 [START_REF] Shreve | From prodigious volcanic degassing to caldera subsidence and quiescence at Ambrym (Vanuatu): the influence of regional tectonics[END_REF].

According to [START_REF] Eissen | Chronology of the historic volcanic activity of the New Hebrides island arc from 1595 to 1991[END_REF], 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 meters), this would result in ∼4 meters of fault slip in a little over 100 years. Assuming the time frame estimated by [START_REF] Robin | Giant tuff cone and 12-km-wide associated caldera at Ambrym Volcano (Vanuatu, New Hebrides Arc)[END_REF] since caldera formation (∼2000 years ago) and stability of Ambrym's activity, we calculate almost 80 meters 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 between (2) intrusion size and fault slip magnitude, would allow for determing 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 non-negligible 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 [START_REF] Geyer | Relationship between caldera collapse and magma chamber withdrawal : An experimental approach[END_REF][START_REF] Geshi | Evaluating volumes for magma chambers and magma withdrawn for caldera collapse[END_REF]. 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 [START_REF] Holohan | Stress evolution during caldera collapse[END_REF].

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 [START_REF] Munro | Caldera morphology in the western Galápagos and implications for volcano eruptive behavior and mechanisms of caldera formation[END_REF]. 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 [START_REF] Amelung | Widespread uplift and 'trapdoor' faulting on Galápagos volcanoes observed with radar interferometry[END_REF][START_REF] Yun | Constraints on magma chamber geometry at Sierra Negra Volcano, Galápagos Islands, based on InSAR observations[END_REF][START_REF] Bell | Caldera resurgence during the 2018 eruption of Sierra Negra volcnao, Galápagos Islands[END_REF]. Much attention has been given to the uplift episodes that cumulate in trap-door faulting along the intra-caldera 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 [START_REF] Jónsson | Stress interaction between magma accumulation and trapdoor faulting on Sierra Negra volcano, Galápagos[END_REF]. Sierra Negra has hosted meter-scale co-eruptive subsidence events in both 2005 and 2018 [START_REF] Geist | The 2005 eruption of Sierra Negra volcano, Galápagos, Ecuador[END_REF][START_REF] Abe | Surface Changes Due to the 2018 Eruption of Sierra Negra Volcano in Galápagos Island Revealed by ALOS-2/PALSAR-2[END_REF][START_REF] Bell | Caldera resurgence during the 2018 eruption of Sierra Negra volcnao, Galápagos Islands[END_REF]. The magnitude of recent co-eruptive subsidence events-5 and 8 meters, respectively-makes them difficult to study with InSAR (S. H. [START_REF] Yun | Interferogram formation in the presence of complex and large deformation[END_REF][START_REF] Casu | Deformation time-series generation in areas characterized by large displacement dynamics: The SAR amplitude pixel-offset SBAS technique[END_REF]. During co-eruptive subsidence events, fault slip has occurred along portions of the sinuous ridge, but there has been no evidence of reactivation of the bounding caldera ring-faults [START_REF] Bell | Caldera resurgence during the 2018 eruption of Sierra Negra volcnao, Galápagos Islands[END_REF].

Future studies concerning subsidence and ring-fault reactivation during moderatesized 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 meters, including subsidence in 2018-2019 and post-eruptive subsidence in 2015, Fig. S14 andS15). 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 pre-eruptive inflation and co-eruptive deflation, which may result in net uplift (i.e., caldera resurgence) or no net displacement over longer periods of time [START_REF] Femina | Deformation of a basaltic shield volcano: Uplift, trapdoor faulting, eruption triggering and subsidence associated with the 2018 eruption of Sierra Negra Volcano, Galápagos[END_REF][START_REF] Galetto | Noneruptive unrest at the caldera of Alcedo volcano (Galápagos Islands) revealed by InSAR data and geodetic modeling[END_REF][START_REF] Bell | Caldera resurgence during the 2018 eruption of Sierra Negra volcnao, Galápagos Islands[END_REF].

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 [START_REF] Shreve | From prodigious volcanic degassing to caldera subsidence and quiescence at Ambrym (Vanuatu): the influence of regional tectonics[END_REF][START_REF] Shreve | Crustal deformation at Ambrym (Vanuatu) imaged with satellite geodesy: constraints on magma storage, migration, and outgassing[END_REF]. 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 modelled from geodetic observations spanning 2004-2020, Fig. S15). The lack of pre-eruptive 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 [START_REF] Coppola | Rheological control on the radiant density of active lava flows and domes[END_REF]. 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 [START_REF] Shreve | Crustal deformation at Ambrym (Vanuatu) imaged with satellite geodesy: constraints on magma storage, migration, and outgassing[END_REF]. 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 inter-collapse uplift, however, more extensive monitoring will be needed to confirm this hypothesis. In addition, Ambrym's caldera formed over thousands of years, a timescale 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 meter 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 Boundary Element Method, 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 Akaike Information Criterion, 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 modelling [START_REF] Hamling | Goldilocks conditions required for earthquakes to trigger basaltic eruptions: Evidence from the 2015 Ambrym eruption[END_REF]. Our study emphasizes the importance of considering mechanical interactions when sources are in close 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 [START_REF] Pascal | On precisely modelling surface deformation due to interacting magma chambers and dykes[END_REF]. Although [START_REF] Pascal | On precisely modelling surface deformation due to interacting magma chambers and dykes[END_REF] 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 -22×10 -10 Pa -1 . The critical fraction of extracted material to reactivate ring-faults at Ambrym has an upper-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 moderatesized eruptions contribute to basaltic caldera development.

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 programme 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. The authors declare no conflict of interests. Data are available on Zenodo with the DOI https://doi.org/10.5281/zenodo.4609496. This is IPGP contribution number 4216. The dike trace represents the upper edge of the quadrangle ("Dtop" in Fig. S7). 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.5×.

Figure 1

 1 Figure 1. a. Tectonic setting of Ambrym volcano, located in the central portion of the New Hebrides subduction zone (pink). AUS is the Australian plate, PAC is the Pacific plate, and DER is the colliding D'Entrecasteaux Ridge. Yellow triangles are active volcanoes, and yellow lines indicate rift zones. The green star indicates the location of the 6.4 Mw earthquake that preceded the 2015 Ambrym eruption. Adapted from Shreve et al. (2019). b. Ambrym island, showing both the caldera and the rift zones. Dates of extra-caldera lava flows associated with rift-zone eruptions are labelled. Adapted from 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 pre-eruptive TanDEM-X DEM and a post-eruptive Pléiades DEM (Fig. S1). Newly opened fissures are denoted by thin dotted lines. e. The cumulative seismic moment and daily SO2 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 Mw earthquake on 19 February 2015 is shown in green. Bounds on the SO2 mass estimate are shown in light blue and are discussed in the Supporting Information Text S1.

  while β m defines the relationship between magmatic pressure change and density change, β m = 1 ρ ∂ρ ∂P . Chamber compressibility β c can be expressed analytically for end-member reservoir shapes (penny-shaped sills (a = b c), spheres (a = b = c), and vertical pipes (a b = c))(Amoruso & Crescentini, 2009), and K. Anderson and Segall (

Figure 2 .

 2 Figure 2. Measurements spanning the February 2015 eruption, including InSAR, MAI, and pixel offsets. The bottom right panel shows the 3D decomposition, zoomed in on the caldera. The arrows represent horizontal displacement and the blue-red colors represent vertical displacement. The arrow in the upper portion of the figure represents 1 meter of eastward displacement. The average and standard deviation of vertical, east-west and north-south displacement within two 1-km wide swath profiles, Profiles A-A' and B-B', are shown in the lower left hand corner of the figure. Positive displacement values indicate uplift, eastward, and northward ground movement. The caldera is outlined with a thick black line, and Marum and Benbow craters are outlined with thin black lines. The red line indicates the fault surface trace.

Figure 3 .

 3 Figure 3. The left hand column shows the ascending and descending InSAR measurements, while the second column shows the modelled synthetic displacement field. The third and fourth columns show the residuals of the modelled and measured displacements, visualized as both unwrapped and wrapped. The data, synthetics, and residuals are shown for a. Inversion 1, which includes a dike, b. Inversion 2, which includes a dike and reservoir, and c. which includes a dike, fault, and reservoir. The bold red traces indicate the sources' projected locations on the map.

Figure 5 .Figure 6 .

 56 Figure 5. The Coulomb stress change on the caldera ring-fault for both Inversion 2 and 3 are shown. An increase in the Coulomb stress indicates that normal dip-slip is promoted. Also shown is the external prestress estimated in the inversion, as well as the final slip distribution on the fault. The total fault slip results from both the external shear stress change and stress transfer from the dike and/or reservoir. In the northern portion of the fault, where the fault-induced surface deformation is measured, the Coulomb stress change is higher in Inversion 3. In Inversion 3, the external stress required on the ring-fault to explain the observed displacement is reduced by almost 40%, indicating that the depressurizing reservoir favors normal slip on the ring-fault.

Table 2 .

 2 A summary of inverted parameters for each inversion. Green cells indicate the parameter is inverted, red cells indicate the parameter is fixed, orange cells are the parameters that are calculated after the inversion, and gray squares indicate the source is not included in that inversion. See FigureS7for an explanation of the listed parameters and TableS5for the explored intervals.

	Parameter	Inversion 1	Inversion 2	Inversion 3
	Dike
	Dip
	Shear
	Bottom Elevation
	Bottom Length
	Twist
	Botang
	D-Top
	Bottom Curvature
	Vertical Length
	Pressure
	Volume Change
	Fault
	Dip
	Shear
	Bottom Elevation
	Bottom Length
	Twist
	Botang
	D-Top
	Bottom Curvature
	Vertical Length
	Shear Stress
	Shear Displacement
	Reservoir
	X0
	Y0
	Semi-major axis
	Semi-minor axis
	Depth
	Strike
	Dip
	Plunge
	Pressure
	Volume Change

Acknowledgments

We thank the Italian Space Agency (ASI Science project No 702) and the Japan Aerospace Exploration Agency (JAXA 6 th 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 was provided through Dinamis Project 2018-147-Sci, ©CNES (2014, 2015, 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 projet ANR-18-CE92-0037. M. Boichu kindly thanks Lieven Clarisse for providing the IASI SO 2 height data and acknowledges the support