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

Tara Shreve, R. Grandin, D. Smittarello, Valérie Cayol, V. Pinel, M. Boichu, Y. Morishita

► To cite this version:

Tara Shreve, R. Grandin, D. Smittarello, Valérie Cayol, V. Pinel, et al. What Triggers Caldera Ring-Fault Subsidence at Ambrym Volcano? Insights From the 2015 Dike Intrusion and Eruption. Journal of Geophysical Research: Solid Earth, 2021, 126 (6), 10.1029/2020JB020277. hal-03353471

HAL Id: hal-03353471

https://uca.hal.science/hal-03353471

Submitted on 24 Sep 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# What triggers caldera ring-fault subsidence at Ambrym volcano? Insights from the 2015 dike intrusion and eruption

T. Shreve<fnref n="1" />*, R. Grandin<fnref n="1" />, D. Smittarello<fnref n="2" />†, V. Cayol<fnref n="3" />, V. Pinel<fnref n="2" />, M. Boichu<fnref n="4" />,
Y. Morishita<fnref n="5" />

## Key Points:

* 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.

<footnote n="1">¹Université de Paris, Institut de physique du globe de Paris, CNRS, F-75005, Paris, France</footnote>

<footnote n="2">²University Grenoble Alpes, University Savoie Mont Blanc, CNRS, IRD, IFSTTAR, ISTerre, Grenoble, France</footnote>

<footnote n="3">³Université Clermont Auvergne, CNRS, IRD, OPGC, Laboratoire Magmas et Volcans, F-6300, Clermont-Ferrand, France</footnote>

<footnote n="4">⁴Université de Lille, UMR 8518 – LOA – Laboratoire d'Optique Atmosphérique, F-59000, Lille, France</footnote>

<footnote n="5">⁵Geospatial Information Authority of Japan, Ibaraki, Japan</footnote>

<footnote>*Currently at Earth and Planets Laboratory, Carnegie Institution for Science, Washington, District of Columbia, USA</footnote>

<footnote>†Currently at the European Center for Geodynamics and Seismology, 19 rue Josy Welter, L-7256 Walferdange, Gd Duchy of Luxembourg</footnote>

<footnote>Corresponding author: Tara Shreve, tshreve@carnegiescience.edu</footnote>
**Abstract**

Surface deformation accompanying dike intrusions is dominated by uplift and horizon-
tal motion directly related to the intrusions. In some cases, it includes subsidence due
to associated magma reservoir deflation. When reservoir deflation is large enough, it can
form, or reactivate pre-existing, caldera ring-faults. Ring-fault reactivation, however, is
rarely observed during moderate-sized eruptions. On February 21st, 2015 at Ambrym
volcano in Vanuatu, a basaltic dike intrusion produced more than 1 meter of co-eruptive
uplift, as measured by InSAR, SAR correlation, and Multiple Aperture Interferometry
(MAI). Here we show that an average of ~40 cm of slip occurred on a normal caldera
ring-fault during this moderate-sized (VEI<3) event, which intruded a volume of ~24×10⁶
m³ and erupted ~9.3×10⁶ m³ of lava (DRE). Using the 3D Mixed Boundary Element
Method, we explore the stress change imposed by the opening dike and the depressur-
izing reservoir on a passive, frictionless fault. Normal fault slip is promoted when stress
is transferred from a depressurizing reservoir beneath one of Ambrym's main craters. Af-
ter estimating magma compressibility, we provide an upper-bound on the critical frac-
tion (f = 7%) of magma extracted from the reservoir to trigger fault slip. We infer that
broad basaltic calderas may form in part by hundreds of subsidence episodes no greater
than a few meters, as a result of magma extraction from the reservoir during moderate-
sized dike intrusions.

**Plain Language Summary**

Many volcanoes feature large depressions, called calderas. Calderas form when enough
magma leaves a deep reservoir, and the solid rock lid above this reservoir can no longer
support its own weight. Caldera faults, or cracks surrounding the reservoir which extend
from the reservoir to Earth's surface, form as the lid collapses. Ambrym volcano (Van-
uatu) has a 12-km wide caldera, and researchers propose it formed during an explosive
eruption 2000 years ago. However, in 2018, Ambrym's caldera sunk along caldera faults
during a non-explosive eruption. This observation questions whether an explosive erup-
tion was necessary to form Ambrym's caldera in the first place. Furthermore, in Febru-
ary 2015, an eruption 10 times smaller than in 2018 also caused the ground to sink along
caldera faults. Utilizing ground motion data obtained from satellite radar systems to model
magma reservoir outflow and fault displacement, we conclude that, in 2015, the ground
sank along caldera faults. This sinking is explained by the removal of as little as <7%
of the stored magma from the reservoir. We therefore propose that Ambrym's wide caldera
may have formed as a result of many frequently occurring medium-sized eruptions. This
challenges the thought that wide calderas mainly form as a result of large eruptions.

**1 Introduction**

Although all calderas host collapse structures, instrumentally recorded caldera ring-
fault activation is rare. A handful of extreme cases of caldera ring-fault activation—catastrophic
caldera collapse—have been observed historically (Neal et al., 2018; Gudmundsson et al.,
2016; Peltier et al., 2009; Geshi et al., 2002). Many of these spectacular events included
noncoherent, or piecemeal, collapse, where upwards propagating faults break up the col-
lapsing piston (e.g., Bárðarbunga (Ágústsdóttir et al., 2019), Kīlauea (Neal et al., 2018;
K. R. Anderson et al., 2019; Segall et al., 2020), Piton de la Fournaise (Peltier et al., 2009)).
In contrast, volcanoes that host shallow and broad magmatic chambers tend to collapse
in a coherent fashion, and the central portion of the collapsing piston stays intact (Roche
et al., 2000). Walker (1984) hypothesized that silicic calderas may form incrementally
by downsagging, which can be accompanied by caldera ring-faulting. The fraction of ma-
terial 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 oc-
cur during frequent eruptions at broad, shallow basaltic calderas. Notably, at Sierra Ne-
gra, 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 ig-
nimbrite deposits (Druitt & Bacon, 1986)) in the literature, despite voluminous histor-
ical eruptions (Reynolds et al., 1995; Munro & Rowland, 1996).

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) (Pinel et al., 2014; Shreve et al., 2019).

Caldera ring-faults can be activated when the reservoir pressure drops beneath litho-
static (i.e., is underpressurized) during magma extraction (Druitt & Sparks, 1984). At
basaltic caldera-rift systems, this extraction occurs primarily through lateral dike intru-
sions (Sigmundsson, 2019). However, in some cases of dike propagation, regardless of whether
or not ring-faults have been activated, surface displacements associated with a depres-
surizing reservoir are not observed geodetically. This may be due to a masking of the
subsidence signal by large displacements related to the dike intrusion (Grandin et al.,
2009). In addition, other factors such as host rock or magma compressibility can reduce
the amplitude of the geodetic signal (Rivalta & Segall, 2008). Accounting for mechan-
ical source interactions allows additional constraints to help identify a depressurizing source.

In February 2015, a dike intrusion fed an eruption at Ambrym volcano and pro-
duced 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.
# 2 Geological Setting

## 2.1 Ambrym Volcano, Vanuatu

Ambrym is a basaltic volcanic island located in the central portion of the New He-
brides 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 (McCall
et al., 1969) (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
(Robin et al., 1993). This view was later challenged by Cronin and Németh (2005), based
on fieldwork reporting the presence of altered mafic deposits around the island, but a
lack of dacites. In December 2018, the first geodetically monitored rift zone intrusion oc-
curred, with >0.4 km³ of magma travelling more than 20 km into the SE rift zone, re-
sulting in a submarine eruption that emitted basaltic pumice onto the nearby shoreline
(Shreve et al., 2019). Coeval with the 2018 rift zone intrusion, the caldera floor subsided
by more than 2 meters. The caldera ring-faults also reactivated, most notably in the north-
ern 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 (Hamling
et al., 2019; Shreve et al., 2019).

This event confirms that the caldera ring-faults of Ambrym are part of an active
fault system, not solely relict structures resulting from a major collapse 2000 years ago
(Robin et al., 1993). Due to these observations of meter-scale caldera subsidence, in con-
junction with historical documentation of reoccurring rift zone intrusions, recent stud-
ies (Hamling et al., 2019; Shreve et al., 2019) have invoked the hypothesis that Ambrym's
caldera developed as a result of hundreds of similar rift zone intrusions, which simulta-
neously 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. How-
ever, caldera faults may reactivate during smaller events. More frequently-occurring moderate-
sized 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 Strom-
bolian explosions within the nested pit craters located in the two main volcanic cones,
Marum and Benbow, near the western caldera rim (see Figure 1c)(Németh & Cronin,
2008). In addition to lava lake activity, occasional intra-caldera fissure eruptions have
occurred, most notably in 1986 and 1988–89 (Eissen et al., 1991) (Fig. 1c). The for-
mer took place in the eastern portion of the caldera, and erupted lava with a slightly el-
evated SiO₂ content (60 wt% SiO₂), compared to eruptive products from the main cones
(~50.5 wt% SiO₂). This suggests various degrees of melt differentiation within the mag-
matic system (Robin et al., 1993; Eissen et al., 1991).

## 2.2 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;
Hamling & Kilgour, 2020). In addition to the main fissure to the SE of Marum, a secondary fissure opened closer to Marum, ~200 m south of Niri Mbwelesu Taten, oriented N133°S, which fed a 900 m-long lava flow (see Figure 1d). This event was preceded by a 6.4 M<sub>w</sub> earthquake on 19 February 2015 at 13:18 UTC with a hypocenter 30 km south-east 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; R. Wright, 2016) and the Middle Infrared Observation of Volcanic Activity (MIROVA) system (Coppola, Laiolo, Cigolini, Delle Donne, & Ripepe, 2016). 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) (Stein et al., 2015; Crawford et al., 2016) and information on the SO₂ plume height derived from Infrared Atmospheric Sounding Interferometer (IASI) observations (Clarisse et al., 2014)). 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⁶ m³ 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 (Rupnik et al., 2018) 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⁶ m³, or ΔV<sub>erupted</sub> = ~9.3 (±0.08) × 10⁶ m³ 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 Bagnardi et al. (2016) and Favalli et al. (2010). 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 conserva-
tively estimate ~40 kt of sulfur dioxide (SO₂) emitted during the eruption (Fig. 1e, Sup-
porting Information Text S3). For this SO₂ mass estimation, we used the OMPS TRM
product, which assumes a center of mass altitude (CMA) of ~7.5 km (Li et al., 2017).
This is in agreement with the retrievals of the SO₂ plume height from IASI observations
(see Text S3). An SO₂ mass estimate of ~40 kt is consistent with degassing of a lava vol-
ume corresponding to the ~9.3×10⁶ m³ lava flow DRE (see Supporting Information Text
S3). After 25 February, SO₂ degassing returned to passive background levels (~7 kt/day,
according to Carn et al. (2017)). These estimates of erupted lava volumes and emitted
gas will be compared with volume changes at depth, as constrained by geodesy. This com-
parison will provide constraints on the mass balance between the erupted and intruded
material, and the overall size of the magma reservoir at depth.
<figure><img src="image_2.png" /><figcaption>**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 M<sub>w</sub> 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 SO<sub>2</sub> mass (as measured by the Ozone Mapping and Profiler Suite, OMPS) emitted throughout the duration of the eruption, which is marked by the orange dotted lines. The timing and focal mechanism of the 6.4 M<sub>w</sub> earthquake on 19 February 2015 is shown in green. Bounds on the SO<sub>2</sub> mass estimate are shown in light blue and are discussed in the Supporting Information Text S1.</figcaption></figure>

## 3 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 InSAR 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.
**3.1 Processing Methods**

**3.1.1 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) (Rosen et al., 2012). Filtering and unwrapping are performed with NSBAS modules (Doin et al., 2011; Rosen et al., 2004; Grandin et al., 2012). 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 (Lindsey et al., 2015; Natsuaki et al., 2016). This was due to an issue with ALOS-2's navigation system, which was fixed on 8 February 2015 (Lindsey et al., 2015). 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 (Ortiz et al., 2007; Natsuaki et al., 2016; Kobayashi et al., 2015). This is implemented with the GSISAR software developed by the Geospatial Information Authority of Japan (Kobayashi et al., 2015; Tobita, 2003). A 30-m mesh Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) Global Digital Elevation Model (GDEM) was used to remove the topographic fringes (Tachikawa et al., 2011). The post-interferogram formation steps (multilooking, filtering, unwrapping and geocoding) are performed with NSBAS modules, as described in the Supporting Information Text S4. The descending interferogram is multilooked 1 time in range and 8 times in azimuth, resulting in a final pixel posting of ~30 m×20 m after geocoding.

**3.1.2 Multiple Aperture Interferometry**

In addition to differential InSAR, which measures satellite line of sight (LOS) displacements, we also process a multiple aperture interferogram (MAI) (Bechor & Zebker, 2006; Liang & Fielding, 2017), in order to derive the along-track, or azimuth, displacement using ascending SM3 images spanning 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 (Goldstein & Werner, 1998), unwrapping, and geocoding with a pixel spacing of ~85 m×70 m.
### 3.1.3 Pixel Offset Tracking

Sub-pixel offset tracking measures surface displacement by finding the cross-correlation peak of two SAR image amplitudes (Michel et al., 1999). Pixel offsets are less accurate than InSAR (Michel et al., 1999). However, they are complementary as they measure both LOS and azimuth displacement, as well as provide measurements in areas where phase is not coherent or surface displacements are large.

For CSK descending acquisitions spanning 13 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).

## 3.2 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 (Grandin et al., 2018). Notable characteristics of the 3D displacement field allow us to hypothesize which deformation sources were active during this event.

## 3.3 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 (Shreve et al., 2019).

## 4 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 (Mogi, 1958), as well as shear or tensile dislocations (Okada, 1985). When multiple sources are involved, the displacement fields produced by each of these sources in an elastic, homogeneous medium are added linearly. Any mechanical source interaction, or stress transfer, between the sources is implicitly neglected.

Given the hypothesis that multiple sources are at play during the 2015 Ambrym eruption, and the geological complexity of caldera systems (pre-existing weaknesses and faults due to caldera collapse (Acocella, 2007)), we choose to proceed with the 3D Mixed Boundary Element Method (BEM) (see Section 4.2 for more details) (Cayol & Cornet, 1997, 1998). This numerical method calculates the surface displacement due to pressure changes within volume sources (or “massive boundaries”) or within fractures with complex geometries (such as curvature, varying dip, etc.), as well as the mechanical source interaction between all 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).

### 4.1 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 (Walstead, 1999; Jónsson et al., 2002), as described in the Supporting Information Text S5, Table S4 and Fig. S5.

Displacement values at subsampled points may be spatially correlated due to phase contributions from atmospheric, ionospheric, or other noise sources (Sudhaus & Jónsson, 2009). To account for this spatially-correlated noise in the interferometric phase, the correlation distance and data variance are used to populate a data covariance matrix, according to $C_d(k_1, k_2) = \sigma_d^2 \exp(-\frac{\|k_1, k_2\|}{a_c})$, where $\sigma_d^2$ is the variance, $a_c$ is the correlation length, and $\|k_1, k_2\|$ is the Euclidian distance between two subsampled pixels $k_1 = (x_1, y_1)$, and $k_2 = (x_2, y_2)$ (Tarantola, 1987). Following Fukushima et al. (2005); Smittarello et al. (2019), we estimated $\sigma_d^2 = 3 \times 10^{-4} \text{ m}^2$ and $a_c = 330 \text{ m}$. The data covariance matrix weights the data in the inversion, taking into account the correlated noise between two pixels.

## 4.2 3D Mixed Boundary Element Method

As previously mentioned, we estimate the geometry and stress change of multiple pressure sources with a 3D Mixed BEM numerical approach (Cayol & Cornet, 1997, 1998). In the 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. Wessel, 2016). Stress changes can either be imposed on the meshed triangular elements of the fractures or massive boundaries, or result from stress transfer from nearby sources, or both.

The Mixed BEM combines the displacement discontinuity and direct displacement methods (Lachat & Watson, 1976; Sokolnikoff, 1956; Crouch, 1976). A description of these methods can be found in the Supporting Information Text S6. Nine parameters control the fracture geometry (Fig. S7) (Fukushima et al., 2010). The fractures are connected to the topography mesh by user-defined surface fissures. In this case, two sets of fissures are organized in an “en echelon” geometry. While their location has a degree of uncertainty, we use high-resolution optical satellite imagery from the Pléiades constellation to map the fissures that opened during the eruption (see Figure 1d). For the caldera ring-fault, the surface trace is determined using the SAR pixel offsets, in conjunction with identifying the incoherent regions in the InSAR measurements.

Displacements and stresses are discretized and interpolated on the boundaries. Thus, the more refined the mesh, the more precise the solutions. There is a 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 (Cayol & Cornet, 1998) and laboratory measurements (Heap et al., 2019). The influence of the chosen shear modulus on the estimated pressure change and reservoir volume will be discussed further in Section 6.2.3.

## 4.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 (Sambridge, 1999b; Fukushima et al., 2005). 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 (Sambridge, 1998; Tridon et al., 2016). Each subsequent iteration picks $N_2$ combinations of parameters, and the $N_2$ computations can be performed in parallel (Fukushima et al., 2005; Smittarello et al., 2019). We take $N_2 = 48$ to optimize this parallelization, which was performed on a cluster using 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 $\chi^2 = (u_o - u_m)^T C_d^{-1} (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, $\sigma_{NLST}$, is reached, based on the standard deviation of the misfit values from the last NLST forward calculations. We set $NLST = N_2$. To calculate posterior probability density functions (PPD's), marginal PPD's, mean model parameters, and model uncertainties, we follow the bayesian inference framework (Tarantola, 1987) described in Sambridge (1999a) and implemented by Fukushima et al. (2005) (Supporting Information Text S7).

The inversion procedure inverts for a phase offset, but does not invert for ramp parameters in the InSAR data that may be due to residual orbital errors (Zebker et al., 1994). Therefore, we mask the displacement signal, fit a ramp to the non-deforming areas, and subtract the ramp from the ascending interferogram using General Mapping Tools (P. Wessel et al., 2019). For the descending interferogram, nearly the entire coherent region is deforming, so we do not estimate and remove a ramp.

## 4.4 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 + \chi^2 + \log|C_d| + N \times \log(2\pi)$, where k is the number of inverted parameters, $\chi^2$ is the misfit, $|C_d|$ is the determinant of the data covariance matrix, and N is the number of data points (Akaike, 1974). In each inversion, we use the same subsampled datasets and data covariance matrix from Section 4.1. Therefore, to compare the AIC of inversions, we can simply compare $AIC = 2k + \chi^2$.

# 5 Inversion Results

## 5.1 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 best-fit model dike is meshed with smaller elements. The final synthetic surface displacements are scaled by a constant to account for discrepancies between computations performed with the coarse and fine meshes (Fukushima et al., 2005). This scalar is ~1 for all inversions (See Table S6, indicating no significant difference in the surface displacements calculated using the coarse and fine meshes). The mean model and 95% uncertainty intervals, as well as the marginal PPD's, are listed in the Supporting Information for all inversions (Table S5 and Figures S9, S10, and S11). The best-fit model dike has a scaled pressure change of $\Delta P_{dike} = 2.2$ MPa and the total volume change is $\Delta V_{dike} = 24 \times 10^6$ m³. 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).

## 5.2 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 (Beauducel et al., 2004; Gudmundsson et al., 2016; Liu et al., 2019). We also only include the portion of the fault that we hypothesize was 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 con-
straint to avoid fracture wall interpenetration. This nonnegativity constraint is enforced
in the Mixed BEM through Lagrange multipliers, as described in Cayol et al. (2014) (from
hereon we call this constraint “preventing interpenetration”). As a result, computation
time of a single forward model increases from 30 seconds to 11 min for the source geom-
etry 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 $\Delta V_{dike} = 23 \times 10^6$ m³, 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).

5.3 Dike, Fault and Reservoir (Inversion 3)

There remain long-wavelength features in the Inversion 2 residuals, postulated to
result from reservoir depressurization, with a similar spatial footprint to the caldera floor
subsidence measured after the Ambrym 2018 rift zone intrusion (see (3) in Section 3.3).
Therefore, a final inversion includes a thin (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 Shreve et al. (2019). Inverting for these parameters would result in an unreasonably
long inverison. As shown by Sambridge (1998), the number of initial computations N₁
needed to adequately explore the parameter space increases exponentially with the num-
ber of inverted parameters. Test inversions also show that the depth cannot be adequately
constrained, and we therefore rely on independent datasets (e.g., Shreve et al. (2019);
Hamling et al. (2019); Hamling and Kilgour (2020); Allard et al. (2015); Legrand et al.
(2005) 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 non-
linear parameters. We approximate the reservoir as a thin oblate spheroid as opposed
to a horizontal fracture because the nonnegativity constraint restricts closing on all frac-
tures 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 ΔVres, reser-
voir volume V and change in pressure ΔPres are inversely proportional (Amoruso & Cres-
centini, 2009), resulting in the high ΔPres obtained by the inversion given the thin reser-
voir. 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 hori-
zontal 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. Pres-
sure change (scaled) of the best-fit model dike is 2.1 MPa, the total volume change of
the dike is ΔV<sub>dike</sub> = 24 × 10<sup>6</sup> m<sup>3</sup>, the total volume change of the reservoir is ΔV<sub>res</sub> =
<table><caption>**Table 1.** The number of free parameters, number of total forward models, AIC (calculated as $AIC = 2k + \chi^2$, where k is the number of inverted parameters and $\chi^2$ is the misfit). Each inversion ended when the standard deviation criterion ($<0.5$) was reached.</caption><thead><tr><th>Deformation sources</th><th>Number of free<br>parameters inverted<br>non-linearly</th><th>Number of total<br>forward models<br>(N<sub>1</sub>/N<sub>2</sub>/N<sub>IT</sub>)</th><th>AIC</th></tr></thead><tbody><tr><td>Dike</td><td>10</td><td>1440/48/463</td><td>18024</td></tr><tr><td>Dike and fault</td><td>13</td><td>9600/48/1789</td><td>14920</td></tr><tr><td>Dike, fault, and reservoir</td><td>15</td><td>9600/48/2819</td><td>10934</td></tr></tbody></table>

$-15 \times 10^6 \text{ 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).

## 5.4 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 (Du et al., 1992). 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 (Shreve et al., 2019; Hamling et al., 2019), as well as depths of Ambrym's shallow reservoir estimated with melt inclusion entrapment pressures and very long period tremors (Allard et al., 2015; Legrand et al., 2005), 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 (Holohan et al., 2011), 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 \times 10^6 \text{ 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.
<table><caption>**Table 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 Figure S7 for an explanation of the listed parameters and Table S5 for the explored intervals.</caption><thead><tr><th>Parameter</th><th>Inversion 1</th><th>Inversion 2</th><th>Inversion 3</th></tr></thead><tbody><tr><td><em>Dike</em></td><td></td><td></td><td></td></tr><tr><td>Dip</td><td></td><td></td><td></td></tr><tr><td>Shear</td><td></td><td></td><td></td></tr><tr><td>Bottom Elevation</td><td></td><td></td><td></td></tr><tr><td>Bottom Length</td><td></td><td></td><td></td></tr><tr><td>Twist</td><td></td><td></td><td></td></tr><tr><td>Botang</td><td></td><td></td><td></td></tr><tr><td>D-Top</td><td></td><td></td><td></td></tr><tr><td>Bottom Curvature</td><td></td><td></td><td></td></tr><tr><td>Vertical Length</td><td></td><td></td><td></td></tr><tr><td>Pressure</td><td></td><td></td><td></td></tr><tr><td>Volume Change</td><td></td><td></td><td></td></tr><tr><td><em>Fault</em></td><td></td><td></td><td></td></tr><tr><td>Dip</td><td></td><td></td><td></td></tr><tr><td>Shear</td><td></td><td></td><td></td></tr><tr><td>Bottom Elevation</td><td></td><td></td><td></td></tr><tr><td>Bottom Length</td><td></td><td></td><td></td></tr><tr><td>Twist</td><td></td><td></td><td></td></tr><tr><td>Botang</td><td></td><td></td><td></td></tr><tr><td>D-Top</td><td></td><td></td><td></td></tr><tr><td>Bottom Curvature</td><td></td><td></td><td></td></tr><tr><td>Vertical Length</td><td></td><td></td><td></td></tr><tr><td>Shear Stress</td><td></td><td></td><td></td></tr><tr><td>Shear Displacement</td><td></td><td></td><td></td></tr><tr><td><em>Reservoir</em></td><td></td><td></td><td></td></tr><tr><td>X0</td><td></td><td></td><td></td></tr><tr><td>Y0</td><td></td><td></td><td></td></tr><tr><td>Semi-major axis</td><td></td><td></td><td></td></tr><tr><td>Semi-minor axis</td><td></td><td></td><td></td></tr><tr><td>Depth</td><td></td><td></td><td></td></tr><tr><td>Strike</td><td></td><td></td><td></td></tr><tr><td>Dip</td><td></td><td></td><td></td></tr><tr><td>Plunge</td><td></td><td></td><td></td></tr><tr><td>Pressure</td><td></td><td></td><td></td></tr><tr><td>Volume Change</td><td></td><td></td><td></td></tr></tbody></table>
# 6 Discussion

## 6.1 Reservoir Depressurization—Not Diking—Induces Faulting

In Section 5.3, we note that Inversion 3 has the lowest AIC, indicating that it is the most likely source configuration. We conclude that contributions from three sources— a dike intrusion, caldera ring-fault, and reservoir—best fit the data, even after the AIC takes into consideration the additional inverted parameters. We will now discuss whether the configuration of three deformation sources is more mechanically consistent than only a dike and ring-fault, providing a supplementary argument, in addition to the AIC, for the presence of a reservoir.

In Inversion 2, we invert for external shear stress change on the caldera ring-fault in order to fit the surface displacements. We consider this an “external” shear stress change because an external stress perturbation must generate the shear stress change, which occurred either prior to the eruption (i.e., the fault was prestressed), during the event itself, or some combination of these events. Static stress transfer from the dike opening during the eruption is the most obvious source of perturbation. 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 (King et al., 1994; Lin & Stein, 2004). The Coulomb stress change is calculated by $\tau_{CSC} = \tau + \mu'\sigma$, where $\tau$ is the shear stress change, $\sigma$ is the normal stress change (positive is unclamping), and $\mu'$ is the effective friction coefficient, where $\mu' = \mu(1 - B)$. B is Skemptons coefficient, which relates pore pressure to confining stress, and we set $\mu' = 0.4$, approximated from laboratory values (Nostro et al., 1998). A positive $\tau_{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 $\tau_{CSC,i}$. The first row in Figure 5 shows the Coulomb stress change on each element of the fault from Inversion 2. We find that the dike opening inhibits normal dip-slip at the fault’s southern end, where surface displacements are dominated by compression perpendicular to the dike plane (Rubin & Pollard, 1988). Normal dip-slip is also inhibited at depth. The Coulomb stress change calculation from Inversion 2 has a slight increase in the northern portion of the fault, in front of the dike where the stress state is dominated by extension. We then calculate the Coulomb stress change from Inversion 3, when the stress perturbation results from a combination of a pressurizing dike and a depressurizing reservoir. There is a positive Coulomb stress change at the fault’s north end, both at the surface and at depth, which promotes normal fault failure (see Figure 5).

Numerical 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 (Buck et al., 2006; Pinel & Jau-
part, 2000; Corbi et al., 2015). 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 re-
verse fault slip at the surface, inconsistent with measured surface displacements. In ad-
dition, 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 pro-
moted in regions of the western rim where we do not observe surface displacement as-
sociated 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 crit-
ical 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. Crit-
ical 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 ex-
ternal 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 (Girona et al., 2014). Therefore, the crit-
ical 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 reser-
voir processes unrelated to magma extraction.

6.2 Caldera Ring-fault Reactivation

**6.2.1 The Influence of Caldera Roof Aspect Ratio**

Caldera ring-faults are formed during caldera collapse. In basaltic systems, calderas develop due to lateral magma propagation which drains a central plumbing system, activating caldera ring-faults and resulting in gradual collapse of the caldera (Sigmundsson, 2019). The plumbing system's geometry influences the critical fraction of magma needed to trigger collapse (defined as $f_{crit} = \frac{-\Delta q_f}{V}$, where $\Delta 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 = \frac{h}{2r}$, where $h$ is reservoir depth and $r$ is the reservoir radius) (Roche et al., 2000; Geshi et al., 2014; Holohan et al., 2011). 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. Anderson et al. (2019) and references within). Ambrym has an especially low roof aspect ratio (<0.4) (Shreve et al., 2019; Hamling et al., 2019), estimated by comparing the depth of the main depressurizing deformation source in 2018 and the width of the caldera surface expression. The 2015 caldera ring-faulting event at Ambrym, although not a full-scale caldera collapse event, allows us to draw bounds on the critical volume fraction needed to reactivate pre-existing caldera ring-faults at a broad, shallow basaltic caldera.

**6.2.2 Pressure Change and Minimum Reservoir Volume**

Using the Mixed BEM and independent obsevations, we adjust the reservoir dimen-
sions to establish bounds on $\Delta P_{res}$ and $V$, after fixing the shear modulus $\mu$ (here we have
$\mu = 2$ GPa, see Section 6.2.3 for a discussion of uncertainties). These bounds will al-
low 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$. Be-
cause $\Delta V_{res} \propto \frac{\Delta P_{res}V}{\mu}$, a tradeoff exists between $\Delta P_{res}$ and $V$ (the full relationship for
various reservoir geometries can be found in Amoruso and Crescentini (2009)). Here, $c$
represents the “height” or “thickness” of the reservoir, which we fix in the inversion be-
cause 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 $\Delta P_{res}$ to obtain a $\Delta V_{res} \approx -15 \times 10^6 \text{ m}^3$ (as estimated in Inver-
sion 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 con-
clude 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 previ-
ous episodes of caldera floor subsidence, such as the subsidence measured in 2015–2017
(Shreve, 2020) and after Ambrym’s 2018 rift zone intrusion (Shreve et al., 2019).

We can put a lower bound on the reservoir’s negative pressure change, because we
do not observe a full-scale caldera collapse creating new ring-faults. A large underpres-
sure Δ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 (Touloukian et al., 1989; Schön,
1996), 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 reser-
voir overpressure before the eruption was at most 15 MPa, given reasonable rock ten-
sile strengths (Sparks, 1997; K. Anderson & Segall, 2013; Delgado et al., 2019). 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³. Given this range, ΔV_res
varies from -13 to -15 × 10⁶ m³.

Although the estimated *V* ranges over multiple orders of magnitude, we note that Allard et al. (2015) estimates Ambrym’s *minimum* reservoir volume to be 0.5 km³. 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₂ flux. To find this minimum, Allard et al. (2015) 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 Carn et al. (2017) (7 kt day⁻¹) and neglect magma recycling into a deeper magmatic system. In this case, the entire magmatic system would contribute to degassing over a period of ~10 years, the time period of the Carn et al. (2017) study. Neglecting effects of magma compressibility, we estimate that 6 km³ of magma was degassed, given ~25,550 kt of SO₂ 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₂ flux span a range of overlapping values <10 km³ (0.02–9.5, 0.5, and
6 km³, respectively). In addition, estimates for the reservoir volume change during the
2018 rift zone intrusion lie between 0.3 and 0.7 km³ (Shreve et al., 2019; Hamling et al.,
2019). If we assume that the reservoir volume change was 0.7 km³ 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³, yet upper bounds on the reservoir size cannot be robustly
constrained in this study.
### 6.2.3 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 (Rivalta & Segall, 2008; Rivalta, 2010). Magma compressibility is one factor that contributes to the discrepancy between estimates of reservoir volume change and the amount of magma extracted from the reservoir (Johnson et al., 2000). Rivalta and Segall (2008); Rivalta (2010) derived the generic formula for $r_v$ based on chamber compressibility $\beta_c$ and magma compressibility $\beta_m$, such that

$$r_v = \frac{\Delta V_s}{\Delta V_c} = 1 + \frac{\beta_m}{\beta_c}, \qquad (1)$$

where $\Delta V_s$ and $\Delta V_c$ are the volume change in the sink (dike + erupted volume) and reservoir, respectively. $\beta_c$ expresses how the host rock responds elastically to pressure change, $\beta_c = \frac{1}{V} \frac{\partial V}{\partial P}$, while $\beta_m$ defines the relationship between magmatic pressure change and density change, $\beta_m = \frac{1}{\rho} \frac{\partial \rho}{\partial P}$. Chamber compressibility $\beta_c$ can be expressed analytically for end-member reservoir shapes (penny-shaped sills ($a = b \gg c$), spheres ($a = b = c$), and vertical pipes ($a \gg b = c$))(Amoruso & Crescentini, 2009), and K. Anderson and Segall (2011) interpolated these estimations using finite element calculations for intermediate aspect ratios.

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 Amoruso and Crescentini (2009), finite element calculations (K. Anderson & Segall, 2011), and the 3D Boundary Element Method. When $a = b \gg 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

$$\beta_c = \frac{3}{\mu} \left( \frac{a}{c} \frac{1}{2\pi} - \frac{1}{5} \right). \qquad (2)$$

Although Equation 2 fails for shallow penny-shaped sills (K. Anderson & Segall, 2011), the “medium” depth (4–5 km b.s.l) of Ambrym’s reservoir allows us to reasonably approximate $\beta_c$ using this expression for lower values of the aspect ratio $c/a$. For the upper-bound, chamber compressibility can be approximated by a sphere within 25% error when aspect ratios $c/a$ range between ~0.6–10 (K. Anderson & Segall, 2011), such that

$$\beta_c = \frac{3}{4\mu}. \qquad (3)$$

Using Equation 2 and 3, it follows that $\beta_c$ ranges from $3.8 \times 10^{-10} – 22 \times 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 \approx 130$ and $V \approx 1$ km$^3$, as concluded in Section 6.2.2). We confirm these analytical approximations by calculating $\beta_c$ using the 3D Boundary Element Method. $\beta_c$ in this case ranges from $3.7 \times 10^{-10} – 18 \times 10^{-10}$. The differences between the analytical and the computational values are consistent with K. Anderson and Segall (2011), and we will proceed with the discussion using the latter range of values.

From Equation 1, given $\Delta V_{dike}$ and $\Delta V_{res}$ estimated in Inversion 3 and $\Delta V_{erupted}$ from Section 2.2, then $r_v = \frac{(\Delta V_{dike} + \Delta V_{erupted})}{\Delta V_{res}} \approx 2.2$ (implying $\beta_m \approx \beta_c$). It follows
that magma compressibility ranges from $4.5 \times 10^{-10} - 22 \times 10^{-10}$ Pa$^{-1}$. Although spanning almost one order of magnitude, we can conclude that this value is higher than that for degassed basalts (Spera, 2000). This estimate is consistent with recent studies of the Ambrym 2015 eruption, which hypothesize that Ambrym's shallow reservoir was pressurized due to bubble nucleation and growth (Hamling & Kilgour, 2020). The existence of gas bubbles within the chamber may explain the relatively high magma compressibility.

Estimates of magma compressibility at Kīlauea during the 2018 lava lake drainage and rift zone eruption ranged from $2 \times 10^{-10} - 15 \times 10^{-10}$ Pa$^{-1}$, also higher than values estimated assuming gas-poor basalt (K. R. Anderson et al., 2019; Segall et al., 2020). K. R. Anderson et al. (2019); Segall et al. (2020) conclude that this may indicate the presence of bubbles in the reservoir, even though the value of $\beta_m$ relies strongly on prior estimates of the shear modulus $\mu$ and $R = \frac{\beta_m}{\beta_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 (Moore et al., 2019). The observations listed above regarding high compressibility magma at Ambrym, Kīlauea, and Erta Ale need to be investigated in more detail before further conclusions can be drawn regarding magma compressibility at volcanoes hosting active lava lakes.

The magma compressibility estimate at Ambrym implies that the reservoir feeding the eruption was not fully degassed. Magma may instead be degassed at shallow levels through persistent lava lake activity. This is consistent with Allard et al. (2015), who concludes that Ambrym's degassing occurs in a "closed system" (Edmonds, 2008). "Closed system" degassing at Ambrym is consistent with the mass balance between the degassing during the eruption and the volume of erupted material. If the melt degasses at very shallow levels, we could assume that there would be no significant degassing through the dike. When comparing the degassed SO$_2$ mass, conservatively 40 kt, which is equivalent to a degassing of $9.5 \times 10^6$ m$^3$ of lava, we indeed notice this is the same order of magnitude as the volume of erupted lava (see Section 2.2). Therefore, we conclude that none of the degassing during this eruption was related to degassing from deeper within the magmatic system, whether from the dike intrusion or the reservoir. The same conclusion was drawn during the December 2018 eruption (Shreve et al., 2019). This balance between the degassed sulfur mass and erupted lava volume has already been documented at other effusive volcanic eruptions (Barnie et al., 2016), yet this is not necessarily the case in explosive eruptions, as emphasized by the "excess" sulfur problem (Wallace, 2001; Kilbride et al., 2016).

In the above discussion, we assume a shear modulus of 2 GPa, a value that is consistent with recent laboratory studies estimating upscaled values of elastic moduli for volcanic rocks (Heap et al., 2019). However, the host rock elastic moduli of a specific region are affected by factors such as porosity, temperature, the presence of microcracks, or confining pressure (Heap et al., 2019). These factors, among others, may cause the shear modulus to vary over more than an order of magnitude. Therefore, we acknowledge that, with no data to constrain the crustal properties at Ambrym, the shear 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 \times 10^{-11}$ Pa$^{-1}$ or as high as $1.8 \times 10^{-8}$ Pa$^{-1}$. Magma compressibility would then vary over a similar range. Comparisons with other magma compressibility estimates should be applied with caution, until studies at Ambrym can further constrain the shear modulus (such as with experimental rock strength estimates using local samples (Heap et al., 2019), or with seismic wave velocities (Ellis et al., 2007; Grandin et al., 2010; Masterlark et al., 2012; Wauthier et al., 2015; Albino et al., 2018)).

An alternative explanation for the discrepancy between the $\Delta V_s$ and $\Delta 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 (Shreve et al., 2019). In addition, the caldera floor subsided in the months following the 2015 eruption, as measured with both ALOS-2 and Sentinel-1 (Fig. S15 and S16). Therefore, at this point, no evidence for magma replenishment can be found.

### 6.2.4 Critical Volume Fractions for Ring-fault Reactivation

Given the estimated $\beta_c$, $\beta_m$, $V$, $\Delta P$, and $\Delta V$, we determine the critical volume fraction needed to trigger caldera ring faulting. Only the western portion of the caldera ring-fault is reactivated, because the reservoir drained during the 2015 event is located beneath the active vents and lava lakes, and this is thus where the largest Coulomb stress change occurred on the ring-fault (See Figure S13b). In the 2018 rift zone intrusion and reservoir drainage event, the northern portion of the caldera ring-fault was reactivated, possibly due to withdrawal of a larger quantity of magma from a more central storage area within Ambrym's magmatic system (Shreve et al., 2019; Hamling et al., 2019).

After doubling the volume change estimates from Inversion 3, as explained in Section 6.1, we assume a volume of reservoir change $-\Delta V_{res,f} = 30 \times 10^6$ m$^3$, a volume of extracted material $-\Delta q = r_v \Delta V_{res,f} = 67 \times 10^6$ m$^3$, and a minimum reservoir volume of $1000 \times 10^6$ m$^3$. Because the reservoir deflation and caldera ring-fault reactivation occurred before the end of the eruption (Fig. S3), the estimated $\Delta 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} = \frac{-\Delta V}{V}$ is thus at most 3%, and, after taking into account the compressibility of the system, the critical fraction $f = \frac{-\Delta q}{V}$ is at most $\sim 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. Geyer et al. (2006) estimate that $f_{exp}$ may be as low as 2% for calderas with the same roof aspect ratio as Ambrym ($r \sim 0.4$). The estimates for $V_{crit}$ and $f$, however, are upper-bounds, considering $\Delta q$ and $\Delta 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).
## 6.3 Implications for Basaltic Caldera Formation

The 2015 eruption at Ambrym is an example of caldera ring-fault reactivation during a moderate-sized eruption. It is possible that localized caldera subsidence may occur during moderate-sized eruptions, which occur on decadal time scales at Ambrym. Caldera-wide subsidence events, on the other hand, occur during rift zone intrusions, because they typically involve a greater volume of magma withdrawal (Shreve et al., 2019). According to Eissen et al. (1991), more than 10 rift zone eruptions have occurred at Ambrym in the past 120 years. Assuming fault slip similar to that measured during the 2018 rift zone intrusion (~0.4 meters), this would result in ~4 meters of fault slip in a little over 100 years. Assuming the time frame estimated by Robin et al. (1993) 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 determining whether incremental caldera subsidence is dominated by larger, rarer events, or by numerous, repeated, moderate events.

The cumulative subsided volume due to these discrete, meter-scale events, which are encouraged due to Ambrym's small roof-aspect ratio (<0.4), is 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 (Geyer et al., 2006; Geshi et al., 2014). For a caldera with a given diameter, a lower roof aspect ratio (shallower reservoir) translates to a ring-fault with a smaller surface area. This lowers the reservoir depressurization needed for caldera ring-fault reactivation. We note that the estimated depressurization and critical volume fraction depend on the physical properties of the fault, such as the friction coefficient, which are not discussed in detail in this study (Holohan et al., 2015).

It is possible that caldera development by ring-faulting during moderate-sized eruptions occurs at other broad and shallow calderas, such as at Sierra Negra in the Galápagos (Munro & Rowland, 1996). Previous studies have concluded that Sierra Negra's reservoir may be a thin ellipsoid 0.4 km tall and 7 km wide at ~2–3 km depth, resulting in a roof aspect ratio of ~0.25–0.4, similar to Ambrym (Amelung et al., 2000; S. Yun et al., 2006; Bell et al., 2021). 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 (Jónsson, 2009). Sierra Negra has hosted meter-scale co-eruptive subsidence events in both 2005 and 2018 (Geist et al., 2008; Abe et al., 2019; Bell et al., 2021). The magnitude of recent co-eruptive subsidence events—5 and 8 meters, respectively—makes them difficult to study with InSAR (S. H. Yun et al., 2007; Casu et al., 2011). 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 (Bell et al., 2021).

Future studies concerning subsidence and ring-fault reactivation during moderate-
sized eruptions at broad mafic calderas, for example, in the Galápagos, Hawaii, or Ice-
land 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 cu-
mulative displacement over geologically significant time scales is net negative or posi-
tive. 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 and S15). In addition, geodetic mea-
surements 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 (La Femina et al., 2019; Galetto et al., 2019;
Bell et al., 2021).

Small episodic collapses may contribute to long-term geomorphological caldera de-
velopment if the amount of inelastic subsidence produced by caldera ring-fault reacti-
vation is greater than any inelastic uplift caused by normal faulting produced during inter-
collapse periods (due to magma replenishment, volatile exsolution, etc.). Ground sub-
sidence at Ambrym may dominate, over decadal time scales, as shown in the uplift/subsidence
ratio measured in the years prior to the December 2018 rift zone intrusion (Shreve et al.,
2019; Shreve, 2020). 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 preced-
ing this event (see cumulative reservoir volume change modelled from geodetic observa-
tions 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 signif-
icant pressurization of the reservoir. Magma replenishment at depth may have occurred
in the years before the rift intrusion, as noted by an increase in thermal output (Coppola
et al., 2013). However, because the system was open, replenishment may have produced
the observed formation of new lava lakes, or the increase in convection vigor and degassing,
as opposed to significant reservoir pressurization producing measurable ground displace-
ment (Shreve, 2020). It is possible that Ambrym is an exceptional case where cumula-
tive 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 hypoth-
esis 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.

# 7 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 (Hamling & Kilgour, 2020). 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 (Pascal et al., 2013). Although Pascal et al. (2013) only discusses the interaction between magma reservoirs and dikes, we have shown that source interactions should also be considered for other geological features, such as faults. In addition to considering source interactions, we also compare the volume change in the dike and the reservoir, estimating a magma compressibility at Ambrym that ranges from $4.5 \times 10^{-10} - 22 \times 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 \text{ km}^3$. However, we acknowledge that a number of parameters influence $f$, including uncertainties on

1. the exact timing of onset of ring-fault slip,

2. the amount of stress change on the ring-fault from prior eruptions, reservoir processes, and/or earthquakes,

3. the frictional parameters of the fault,

4. the geometry of the system, especially the connection between the fault and the reservoir, and

5. the reservoir volume's unknown upper bound.

Continued geodetic monitoring of Ambrym, in conjunction with widespread observations at basaltic calderas around the world, will allow us to further constrain how moderate-sized eruptions contribute to basaltic caldera development.

## Acknowledgments

We thank the Italian Space Agency (ASI Science project No 702) and the Japan Aerospace Exploration Agency (JAXA 6th Research Agreement No 3245) for providing access to the radar imagery used in this study. A portion of the ALOS-2 data was provided under a cooperative research contract between Geospatial Information Authority of Japan and JAXA. The ownership of ALOS-2 data belongs to JAXA. The Pléiades data 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
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.

## References

Abe, T., Ohki, M., & Tadono, T. (2019). Surface Changes Due to the 2018 Eruption of Sierra Negra Volcano in Galápagos Island Revealed by ALOS-2/PALSAR-2. *IGARSS 2019 - 2019 IEEE International Geoscience and Remote Sensing Symposium*, 9334-9337. doi: 10.1109/igarss.2019.8900424

Acocella, V. (2007). Understanding caldera structure and development: An overview of analogue models compared to natural calderas. *Earth-Science Reviews*, 85(3-4), 125-160. doi: 10.1016/j.earscirev.2007.08.004

Ágústsdóttir, T., Winder, T., Woods, J., White, R. S., Greenfield, T., & Brandsdóttir, B. (2019). Intense seismicity during the 2014-2015 Bárðarbunga-Holuhraun rifting event, Iceland, reveals the nature of dike-induced earthquakes and caldera collapse mechanisms. *Journal of Geophysical Research: Solid Earth*, 124(8), 8331-8357. doi: 10.1029/2018JB016010

Akaike, H. (1974). A new look at the statistical model identification. *IEEE Transactions on Automatic Control*, 19(6), 716-723. doi: 10.1109/TAC.1974.1100705

Albino, F., Amelung, F., & Gregg, P. (2018). The Role of Pore Fluid Pressure on the Failure of Magma Reservoirs: Insights From Indonesian and Aleutian Arc Volcanoes. *Journal of Geophysical Research: Solid Earth*, 123(2), 1328-1349. doi: 10.1002/2017JB014523

Allard, P., Aiuppa, A., Bani, P., Métrich, N., Bertagnini, A., Gauthier, P., ... Garaebiti, E. (2015). Prodigious emission rates and magma degassing budget of major, trace and radioactive volatile species from Ambrym basaltic volcano, Vanuatu island Arc. *Journal of Volcanology and Geothermal Research*, 304, 378-402. doi: 10.1016/j.jvolgeores.2015.08.022

Amelung, F., Zebker, H., & Segall, P. (2000). Widespread uplift and ‘trapdoor’ faulting on Galápagos volcanoes observed with radar interferometry. *October*, 407(October), 993-996. doi: 10.1038/nature01526.1

Amoruso, A., & Crescentini, L. (2009). Shape and volume change of pressurized ellipsoidal cavities from deformation and seismic data. *Journal of Geophysical Research: Solid Earth*, 114(2), 1-8. doi: 10.1029/2008JB005946

Anderson, K., & Segall, P. (2011). Physics-based models of ground deformation and extrusion rate at effusively erupting volcanoes. *Journal of Geophysical Research: Solid Earth*, 116(7), 1-20. doi: 10.1029/2010JB007939

Anderson, K., & Segall, P. (2013). Bayesian inversion of data from effusive volcanic eruptions using physics-based models: Application to Mount St. Helens 2004-
2008. *Journal of Geophysical Research: Solid Earth*, 118(5), 2017–2037. doi: 10.1002/jgrb.50169

Anderson, K. R., Johanson, I. A., Patrick, M. R., Gu, M., Segall, P., Poland, M. P., ... Miklius, A. (2019). Magma reservoir failure and the onset of caldera collapse at Kīlkōlcano in 2018. *Science* (6470). doi: 10.1126/science.aaz1822

Bagnardi, M., González, P. J., & Hooper, A. (2016). High-resolution digital elevation model from tri-stereo Pleiades-1 satellite imagery for lava flow volume estimates at Fogo Volcano. *Geophysical Research Letters*, 43(12), 6267–6275. doi: 10.1002/2016GL069457

Barnie, T. D., Keir, D., Hamling, I., Hofmann, B., Belachew, M., Carn, S., ... Wright, T. (2016). A multidisciplinary study of the final episode of the Manda Hararo dyke sequence, Ethiopia, and implications for trends in volcanism during the rifting cycle. *Geological Society, London, Special Publications*, 420(1), 149–163. doi: 10.1144/SP420.6

Beauducel, F., De Natale, G., Obrizzo, F., & Pingue, F. (2004). 3-D modelling of Campi Flegrei ground deformations: Role of Caldera boundary discontinuities. *Pure and Applied Geophysics*, 161(7), 1329–1344. doi: 10.1007/s00024-004-2507-4

Bechor, N. B., & Zebker, H. A. (2006). Measuring two-dimensional movements using a single InSAR pair. *Geophysical Research Letters*, 33(16), 1–5. doi: 10.1029/2006GL026883

Bell, A., La Femina, P., Ruiz, M., Amelung, F., Bagnardi, M., Bean, C., ... Stock, M. (2021). Caldera resurgence during the 2018 eruption of Sierra Negra volcanao, Galápagos Islands. *Nature Communications*, 12(1397), 1–9. doi: 10.1038/s41467-021-21596-4

Buck, W. R., Einarsson, P., & Brandsdóttir, B. (2006). Tectonic stress and magma chamber size as controls on dike propagation: Constraints from the 1975–1984 Krafla rifting episode. *Journal of Geophysical Research: Solid Earth*, 111(12), 1–15. doi: 10.1029/2005JB003879

Carn, S. A., Fioletov, V. E., Mclinden, C. A., Li, C., & Krotkov, N. A. (2017). A decade of global volcanic SO₂ emissions measured from space. *Scientific Reports*, 7, 1–12. doi: 10.1038/srep44095

Casu, F., Manconi, A., Pepe, A., & Lanari, R. (2011). Deformation time-series generation in areas characterized by large displacement dynamics: The SAR amplitude pixel-offset SBAS technique. *IEEE Transactions on Geoscience and Remote Sensing*, 49(7), 2752–2763. doi: 10.1109/TGRS.2010.2104325

Cayol, V., Catry, T., Michon, L., Chaput, M., Famin, V., Bodart, O., ... Romagnoli, C. (2014). Sheared sheet intrusions as mechanism for lateral flank displacement on basaltic volcanoes: Applications to Réunion Island volcanoes. *Journal of Geophysical Research: Solid Earth*, 119.

Cayol, V., & Cornet, F. H. (1997). 3D mixed boundary elements for elastostatic deformation field analysis. *International journal of rock mechanics and mining*, 34(2), 275–287. doi: 10.1016/S0148-9062(96)00035-6

Cayol, V., & Cornet, F. H. (1998). Three-dimensional modeling of the 1983–1984
eruption at Piton de la Fournaise Volcano, Réunion Island. *Journal of Geophysical Research*, 103, 18,025–18,037.

Clarisse, L., Coheur, P. F., Theys, N., Hurtmans, D., & Clerbaux, C. (2014). The 2011 Nabro eruption, a SO₂ plume height analysis using IASI measurements. *Atmospheric Chemistry and Physics*, 14(6), 3095–3111. doi: 10.5194/acp-14-3095-2014

Coppola, D., Laiolo, M., & Cigolini, C. (2016). Fifteen years of thermal activity at Vanuatu's volcanoes (2000 - 2015) revealed by MIROVA. *Journal of Volcanology and Geothermal Research*, 322, 6–19. doi: 10.1016/j.jvolgeores.2015.11.005

Coppola, D., Laiolo, M., Cigolini, C., Delle Donne, D., & Ripepe, M. (2016). Enhanced volcanic hot-spot detection using MODIS IR data: Results from the MIROVA system. *Geological Society Special Publication*, 426(1), 181–205. doi: 10.1144/SP426.5

Coppola, D., Laiolo, M., Piscopo, D., & Cigolini, C. (2013). Rheological control on the radiant density of active lava flows and domes. *Journal of Volcanology and Geothermal Research*, 249, 39–48. doi: 10.1016/j.jvolgeores.2012.09.005

Corbi, F., Rivalta, E., Pinel, V., Maccaferri, F., Bagnardi, M., & Acocella, V. (2015). How caldera collapse shapes the shallow emplacement and transfer of magma in active volcanoes. *Earth and Planetary Science Letters*, 431, 287–293. doi: 10.1016/j.epsl.2015.09.028

Crawford, A. M., Stunder, B. J., Ngan, F., & Pavolonis, M. J. (2016). Initializing HYSPLIT with satellite observations of volcanic ash: A case study of the 2008 Kasatochi eruption. *Journal of Geophysical Research*, 121(18), 10,786–10,803. doi: 10.1002/2016JD024779

Cronin, S. J., & Németh, K. (2005). Where are the giant tuff cone and ignimbrites of Ambrym? A more conventional story of mafic volcanism at Ambrym Volcano, Vanutu. In *Geological society of new zealand 50th annual conference* (pp. 21–22). Kaikoura, New Zealand: Geological Society of New Zealand.

Crouch, S. L. (1976). Solution of plane elasticity problems by the displacement discontinuity method. *International Journal for Numerical Methods in Engineering*, 10(2), 301–343.

Delgado, F., Kubanek, J., Anderson, K., Lundgren, P., & Pritchard, M. (2019). Physicochemical models of effusive rhyolitic eruptions constrained with InSAR and DEM data: A case study of the 2011–2012 Cordón Caulle eruption. *Earth and Planetary Science Letters*, 524, 115736. doi: 10.1016/j.epsl.2019.115736

Doin, M.-P., Lodge, F., Guilloso, S., Jolivet, R., Lasserre, C., Ducret, G., ... Pinel, V. (2011). Presentation of the small baseline NSBAS processing chain on a case example: the Etna deformation monitoring from 2003 to 2010 using Envisat data. In *Proceedings of the esa 'fringe 2011 workshop', frascati, italy*, (19-23 september 2011) (pp. 19–23).

Druitt, T., & Bacon, C. (1986). Lithic breccia and ignimbrite erupted during the collapse of Crater Lake Caldera, Oregon. *Journal of Volcanology and Geothermal Research*, 29(1-4), 1–32. doi: 10.1016/0377-0273(86)90038-7
Druitt, T., & Sparks, R. (1984). On the formation of calderas during ignimbrite eruptions. *Nature*, **310**(5979), 679–681. doi: 10.1038/310679a0

Du, Y., Aydin, A., & Segall, P. (1992). Comparison of various inversion techniques as applied to the determination of a geophysical deformation model for the 1983 Borah Peak earthquake. *Society, 82*(4), 1840–1866.

Edmonds, M. (2008). New geochemical insights into volcanic degassing. *Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences*, **366**(1885), 4559–4579. doi: 10.1098/rsta.2008.0185

Eissen, J. P., Blot, C., & Louat, R. (1991). *Chronology of the historic volcanic activity of the New Hebrides island arc from 1595 to 1991* (Tech. Rep.). Nouméa: ORSTOM.

Ellis, S., Wilson, C., Bannister, S., Bibby, H., Heise, W., Wallace, L., & Patterson, N. (2007). A future magma inflation event under the rhyolitic Taupo volcano, New Zealand: Numerical models based on constraints from geochemical, geological, and geophysical data. *Journal of Volcanology and Geothermal Research*, **168**(1–4), 1–27. doi: 10.1016/j.jvolgeores.2007.06.004

Favalli, M., Fornaciai, A., Mazzarini, F., Harris, A., Neri, M., Behncke, B., ...
Boschi, E. (2010). Evolution of an active lava flow field using a multitemporal
LIDAR acquisition. *Journal of Geophysical Research: Solid Earth*, **115**(11),
1–17. doi: 10.1029/2010JB007463

Fukushima, Y., Cayol, V., & Durand, P. (2005). Finding realistic dike models from interferometric synthetic aperture radar data: The February 2000 eruption at Piton de la Fournaise. *Journal of Geophysical Research: Solid Earth*, **110**(3), 1–15. doi: 10.1029/2004JB003268

Fukushima, Y., Cayol, V., Durand, P., & Massonnet, D. (2010). Evolution of magma conduits during the 1998-2000 eruptions of Piton de la Fournaise volcano, Réunion Island. *Journal of Geophysical Research: Solid Earth*, **115**(10), 1–21. doi: 10.1029/2009JB007023

Galetto, F., Bagnardi, M., Acocella, V., & Hooper, A. (2019). Nonruptive unrest at the caldera of Alcedo volcano (Galápagos Islands) revealed by InSAR data and geodetic modeling. *Journal of Geophysical Research: Solid Earth*, **124**(4), 3365–3381. doi: 10.1029/2018JB017103

Geist, D. J., Harpp, K. S., Naumann, T. R., Poland, M., Chadwick, W. W., Hall, M., & Rader, E. (2008). The 2005 eruption of Sierra Negra volcano, Galápagos, Ecuador. *Bulletin of Volcanology*, **70**(6), 655–673. doi: 10.1007/s00445-007-0160-3

Geshi, N., Ruch, J., & Acocella, V. (2014). Evaluating volumes for magma chambers and magma withdrawn for caldera collapse. *Earth and Planetary Science Letters*, **396**, 107–115. doi: 10.1016/j.epsl.2014.03.059

Geshi, N., Shimano, T., Chiba, T., & Nakada, S. (2002). Caldera collapse during the 2000 eruption of Miyakejima Volcano, Japan. *Bulletin of Volcanology*, **64**(1), 55–68. doi: 10.1007/s00445-001-0184-z

Geyer, A., Folch, A., & Martí, J. (2006). Relationship between caldera collapse and magma chamber withdrawal: An experimental approach. *Journal of Volcanol-
*ogy and Geothermal Research*, **157**, 375–386. doi: 10.1016/j.jvolgeores.2006.05.001

Girona, T., Costa, F., Newhall, C., & Taisne, B. (2014). On depressurization of volcanic magma reservoirs by passive degassing. *Journal of Geophysical Research: Solid Earth*, **119**, 8667–8687. doi: 10.1002/2014JB011368.Received

Goldstein, R. M., & Werner, C. L. (1998). Radar Interferogram Filtering for Geophysical Applications. *Geophysical Research Letters*, **25**(21), 4035–4038. doi: https://doi.org/10.1029/1998GL900033

Grandin, R., Doin, M.-P., Bollinger, L., Pinel-Puysségur, B., Ducret, G., Jolivet, R., & Sapkota, S. N. (2012). Long-term growth of the Himalaya inferred from interseismic InSAR measurement. *Geology*, **40**(12), 1059. doi: 10.1130/G33154.1

Grandin, R., Lauer, B., Klinger, Y., & Jolivet, R. (2018). 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 (2013, M<sub>w</sub> 7.7, Pakistan). *American Geophysical Union Fall Meeting*. Retrieved from https://agu.confex.com/agu/fm18/meetingapp.cgi/Paper/452841

Grandin, R., Socquet, A., Binet, R., Klinger, Y., Jacques, E., De Chabalier, J. B., ... Pinzuti, P. (2009). September 2005 Manda hararo-dabbahu rifting event, Afar (Ethiopia): Constraints provided by geodetic data. *Journal of Geophysical Research: Solid Earth*, **114**(8). doi: 10.1029/2008JB005843

Grandin, R., Socquet, A., Jacques, E., Mazzoni, N., De Chabalier, J., & King, G. (2010). 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. *Journal of Geophysical Research: Solid Earth*, **115**(10), 2005–2009. doi: 10.1029/2009JB000815

Gudmundsson, M. T., Jónsdóttir, K., Hooper, A., Holohan, E. P., Halldórsson, S. A., Ófeigsson, B. G., ... Aiuppa, A. (2016). Gradual caldera collapse at Bárdarbunga volcano, Iceland, regulated by lateral magma outflow. *Science*, **353**(6296), aaf8988. doi: 10.1126/science.aaf8988.sciencemag.org

Hamling, I. J., Cevuard, S., & Garaebiti, E. (2019). Large-scale drainage of a complex magmatic system: Observations from the 2018 eruption of Ambrym volcano, Vanuatu. *Geophysical Research Letters*, **46**(9), 4609–4617. doi: 10.1029/2019GL082606

Hamling, I. J., & Kilgour, G. (2020). Goldilocks conditions required for earthquakes to trigger basaltic eruptions: Evidence from the 2015 Ambrym eruption. *Science Advances*, **6**(14). doi: 10.1126/sciadv.aaz5261

Heap, M. J., Villeneuve, M., Albino, F., Farquharson, J. I., Brothelande, E., Amelung, F., ... Baud, P. (2019). Towards more realistic values of elastic moduli for volcano modelling. *Journal of Volcanology and Geothermal Research*, 106684. doi: 10.1016/j.jvolgeores.2019.106684

Holohan, E. P., Schöpfer, M. P., & Walsh, J. J. (2011). Mechanical and geometric controls on the structural evolution of pit crater and caldera subsidence. *Journal of Geophysical Research: Solid Earth*, **116**(7). doi: 10.1029/2010JB008032
Holohan, E. P., Schöpfer, M. P. J., & Walsh, J. J. (2015). Stress evolution during caldera collapse. *Earth and Planetary Science Letters*, **421**, 139–151.

Johnson, D. J., Sigmundsson, F., & Delaney, P. T. (2000). 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. *Bulletin of Volcanology* **(7)**, 491–493. doi: 10.1007/s004450050006

Jónsson, S. (2009). Stress interaction between magma accumulation and trapdoor faulting on Sierra Negra volcano, Galápagos. *Tectonophysics*, **471**(1-2), 36-44.
doi: 10.1016/j.tecto.2008.08.005

Jónsson, S., Zebker, H., Segall, P., & Amelung, F. (2002). Fault slip distribution of the 1999 M<sub>w</sub> 7.1 Hector Mine, California, earthquake, estimated from satellite radar and GPS measurements. *Bulletin of the Seismological Society of America*, **92**(4), 1377–1389. doi: 10.1785/0120000922

Kilbride, B. M. C., Edmonds, M., & Biggs, J. (2016). Observing eruptions of gas-rich compressible magmas from space. *Nature Communications*, **7**, 1–8. doi: 10.1038/ncomms13744

King, G. C. P., Stein, R. S., & Lin, J. (1994). Static Stress Change and the Triggering of Earthquakes. *Bulletin of the Geological Society of America* (March), 46.

Kobayashi, T., Morishita, Y., & Yarai, H. (2015). Detailed crustal deformation and fault rupture of the 2015 Gorkha earthquake, Nepal, revealed from ScanSAR-based interferograms of ALOS-2. *Earth, Planets and Space*, **67**(1). doi: 10.1186/s40623-015-0359-z

La Femina, P. C., Bell, A. F., Higgins, M., Meier, N. J., Hernandez, S., Bean, C. J., ... Grannell, J. (2019, Dec). Deformation of a basaltic shield volcano: Uplift, trapdoor faulting, eruption triggering and subsidence associated with the 2018 eruption of Sierra Negra Volcano, Galápagos. In *Agu fall meeting abstracts* (p. V11B-01).

Lachat, J. C., & Watson, J. O. (1976). Effective numerical treatment of boundary integral equations: A formulation for three-dimensional elastostatics. *International Journal for Numerical Methods in Engineering*, **10**(5), 991–1005. doi: 10.1002/nme.1620100503

Legrand, D., Rouland, D., Frogneux, M., Carniel, R., Charley, D., Roult, G., & Robin, C. (2005). Interpretation of very long period tremors at Ambrym volcano, Vanuatu, as quasi-static displacement field related to two distinct magmatic sources. *Geophysical Research Letters*, **32**(6), 1–5. doi: 10.1029/2004GL021968

Li, C., Krotkov, N. A., Carn, S., Zhang, Y., Spurr, R. J., & Joiner, J. (2017). New-generation NASA Aura Ozone Monitoring Instrument (OMI) volcanic SO₂ dataset: Algorithm description, initial results, and continuation with the Suomi-NPP Ozone Mapping and Profiler Suite (OMPS). *Atmospheric Measurement Techniques*, **10**(2), 445–458. doi: 10.5194/amt-10-445-2017

Liang, C., & Fielding, E. J. (2017). Measuring Azimuth Deformation with L-Band ALOS-2 ScanSAR Interferometry. *IEEE Transactions on Geoscience and Re-
mote Sensing, 55(5), 2725-2738. doi: 10.1109/TGRS.2017.2653186

Lin, J., & Stein, R. S. (2004). Stress triggering in thrust and subduction earthquakes and stress interaction between the southern San Andreas and nearby thrust and strike-slip faults. *Journal of Geophysical Research: Solid Earth, 109*(B2), 1-19. doi: 10.1029/2003jb002607

Lindsey, E. O., Natsuaki, R., Xu, X., Shimada, M., Hashimoto, M., Melgar, D., & Sandwell, D. T. (2015). Line-of-sight displacement from ALOS-2 interferometry: $M_w$ 7.8 Gorkha Earthquake and $M_w$ 7.3 aftershock. *Geophysical Research Letters*, 42(16), 6655-6661. doi: 10.1002/2015GL065385

Liu, Y.-K., Ruch, J., Vasyura-Bathke, H., & Jónsson, S. (2019). Influence of ring faulting in localizing surface deformation at subsiding calderas. *Earth and Planetary Science Letters*, 526, 115784. doi: 10.1016/j.epsl.2019.115784

Masterlark, T., Feigl, K., Haney, M., Stone, J., Thurber, C., & Ronchin, E. (2012). Nonlinear estimation of geometric parameters in FEMs of volcano deformation: Integrating tomography models and geodetic data for Okmok volcano, Alaska. *Journal of Geophysical Research: Solid Earth*, 117(2), 1-17. doi: 10.1029/2011JB008811

McCall, G., Lemaitre, R., Malahoff, A., Robinson, G., & Stephenson, P. (1969). The Geology and Geophysics of the Ambrym Caldera, New Hebrides. In *Symposium volcanoes and their roots* (pp. 681-696). Oxford, England.

Michel, R., Avouac, J.-P., & Taboury, J. (1999). Measuring ground displacements from SAR amplitude images: application to the Landers earthquake. *Geophysical Research Letters*, 26(7), 875-878.

Mogi, K. (1958). Relations between the eruptions of various volcanoes and the deformations of the ground surfaces around them. *Bulletin of the Earthquake Research Institute*, 36, 99-134.

Moore, C., Wright, T., Hooper, A., & Biggs, J. (2019). The 2017 Eruption of Erta 'Ale Volcano, Ethiopia: Insights into the Shallow Axial Plumbing System of an Incipient Mid-Ocean Ridge. *Geochemistry, Geophysics, Geosystems*, 20, 5727-5743. doi: 10.1029/2019GC008692

Munro, D. C., & Rowland, S. K. (1996). Caldera morphology in the western Galápagos and implications for volcano eruptive behavior and mechanisms of caldera formation. *Journal of Volcanology and Geothermal Research*, 72, 85-100.

Natsuaki, R., Nagai, H., Motohka, T., Ohki, M., Watanabe, M., B. Thapa, R., ... Suzuki, S. (2016). SAR interferometry using ALOS-2 PALSAR-2 data for the $M_w$ 7.8 Gorkha, Nepal earthquake. *Earth, Planets and Space*, 68(1). doi: 10.1186/s40623-016-0394-4

Neal, C. A., Brantley, S., Antolik, L., Babb, J., Burgess, M., Calles, K., ... Fisher, G. (2018). The 2018 rift eruption and summit collapse of Kīlauea volcano. *Science*, 1-14. doi: 10.1126/science.aav7046

Németh, K., & Cronin, S. J. (2008). Volcanic craters, pit craters and high-level magma-feeding systems of a mafic island-arc volcano: Ambrym, Vanuatu, South Pacific. *Geological Society, London, Special Publications*, 302(1), 87-102. doi: 10.1144/SP302.6
Nostro, C., Stein, S., Cocco, M., Belardinelli, M. E., & Marzocchi, W. (1998). Two-way coupling between Vesuvius eruptions and southern Apennine earthquakes, Italy, by elastic stress transfer. *Journal of Geophysical Research, 103*(B10), 24,487-24,504.

Okada, Y. (1985). Surface Deformation Due to Shear and Tensile Faults in a Half-Space. *Bulletin of the Seismological Society of America, 75*(4), 1135-1154. doi: 10.1016/0148-9062(86)90674-1

Ortiz, A. B., Member, S., & Zebker, H. (2007). ScanSAR-to-Stripmap Mode Interferometry Processing Using ENVISAT / ASAR Data. *IEEE Transactions on Geoscience and Remote Sensing, 45*(11), 3468-3480. doi: 10.1109/TGRS.2007.895970

Pascal, K., Neuberg, J., & Rivalta, E. (2013). On precisely modelling surface deformation due to interacting magma chambers and dykes. *Geophysical Journal International, 196*(1), 253-278. doi: 10.1093/gji/ggt343

Peltier, A., Staudacher, T., Bachèlery, P., & Cayol, V. (2009). Formation of the April 2007 caldera collapse at Piton de La Fournaise volcano: Insights from GPS data. *Journal of Volcanology and Geothermal Research, 184*(1-2), 152-163. doi: 10.1016/j.jvolgeores.2008.09.009

Pinel, V., & Jaupart, C. (2000). The effect of edifice load on magma ascent beneath a volcano. *Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 358*(1770), 1515-1532. doi: 10.1098/rsta.2000.0601

Pinel, V., Poland, M. P., & Hooper, A. (2014). Volcanology: Lessons learned from Synthetic Aperture Radar imagery. *Journal of Volcanology and Geothermal Research, 289*, 81-113. doi: 10.1016/j.jvolgeores.2014.10.010

Reynolds, R. W., Geist, D., & Kurz, M. D. (1995). Physical volcanology and structural development of Sierra Negra volcano, Isabela Island, Galápagos archipelago. *Geological Society of America Bulletin, 107*(12), 1398-1410. doi: 10.1130/0016-7606(1995)107<1398:PVASDO>2.3.CO;2

Rivalta, E. (2010). Evidence that coupling to magma chambers controls the volume history and velocity of laterally propagating intrusions. *Journal of Geophysical Research: Solid Earth, 115*(7), 1-13. doi: 10.1029/2009JB006922

Rivalta, E., & Segall, P. (2008). Magma compressibility and the missing source for some dike intrusions. *Geophysical Research Letters, 35*(4), 1-5. doi: 10.1029/2007GL032521

Robin, C., Eissen, J. P., & Monzier, M. (1993). Giant tuff cone and 12-km-wide associated caldera at Ambrym Volcano (Vanuatu, New Hebrides Arc). *Journal of Volcanology and Geothermal Research, 55*(3-4), 225-238. doi: 10.1016/0377-0273(93)90039-T

Roche, O., & Druitt, T. (2001). Onset of caldera collapse during ignimbrite eruptions. *Earth and Planetary Science Letters, 191*(3-4), 191-202. doi: 10.1016/S0012-821X(01)00428-9

Roche, O., Druitt, T. H., & Merle, O. (2000). Experimental study of caldera formation. *Journal of Geophysical Research, 105*(B1), 395-416.
Roche, O., Van Wyk De Vries, B., & Druitt, T. H. (2001). Sub-surface structures and collapse mechanisms of summit pit craters. *Journal of Volcanology and Geothermal Research*, 105(1-2), 1-18. doi: 10.1016/S0377-0273(00)00248-1

Rosen, P. A., Gurrola, E., Sacco, G. F., & Zebker, H. (2012, April). The InSAR scientific computing environment. In *Eusar 2012; 9th european conference on synthetic aperture radar* (p. 730-733).

Rosen, P. A., Hensley, S., Peltzer, G., & Simons, M. (2004). Updated Repeat Orbit Interferometry Package Released. *Eos*, 85(5), 47.

Rubin, A. M., & Pollard, D. D. (1988). Dike-induced faulting in rift zones of Iceland and Afar. *Geology*, 16, 413-417.

Rupnik, E., Pierrot-Deseilligny, M., & Delorme, A. (2018). 3D reconstruction from multi-view VHR-satellite images in MicMac. *ISPRS Journal of Photogrammetry and Remote Sensing*, 139, 201-211. doi: 10.1016/j.isprsjprs.2018.03.016

Sambridge, M. (1998). Exploring multidimensional landscapes without a map. *Inverse Problems* (14), 427-440.

Sambridge, M. (1999a). Geophysical inversion with a neighbourhood algorithm— II. Appraising the ensemble. *Geophysical Journal International*, 138(2), 727-746. doi: 10.1046/j.1365-246X.1999.00876.x

Sambridge, M. (1999b). Geophysical inversion with a neighbourhood algorithm— I. Searching a parameter space. *Geophysical Journal International*, 138, 479-494. doi: 10.1046/j.1365-246x.1999.00900.x

Schön, J. (1996). Fundamentals and principles of petrophysics, handbook of geo-physical exploration, seismic exploration. In K. Helbig & S. Treitel (Eds.), *Physical properties of rocks*. (p. 583).

Segall, P., Anderson, K., Pulvirenti, F., Wang, T., & Johanson, I. (2020). Caldera collapse geometry revealed by near-field GPS displacements at Kīlauea volcano in 2018. *Geophysical Research Letters* (May 2018), 1-8. doi: http://doi.org/10.1029/2020GL088867

Shreve, T. (2020). *Crustal deformation at Ambrym (Vanuatu) imaged with satellite geodesy: constraints on magma storage, migration, and outgassing* (Unpublished doctoral dissertation). Institut de Physique du Globe de Paris.

Shreve, T., Grandin, R., Boichu, M., Garaebiti, E., Moussallam, Y., Ballu, V., ...
Pelletier, B. (2019). From prodigious volcanic degassing to caldera subsidence and quiescence at Ambrym (Vanuatu): the influence of regional tectonics. *Scientific Reports*, 9(18868). doi: https://doi.org/10.1038/s41598-019-55141-7

Sigmundsson, F. (2019). Calderas collapse as magma flows into rifts. *Science*, 366(6470), 1200-1201. doi: 10.1126/science.aaz7126

Smittarello, D., Cayol, V., Pinel, V., Peltier, A., Froger, J. L., & Ferrazzini, V. (2019). Magma Propagation at Piton de la Fournaise From Joint Inversion of InSAR and GNSS. *Journal of Geophysical Research: Solid Earth*, 124(2), 1361-1387. doi: 10.1029/2018JB016856

Sokolnikoff, I. (1956). *Mathematical theory of elasticity*. New York: McGraw-Hill.

Sparks, R. S. (1997). Causes and consequences of pressurisation in lava dome eruptions. *Earth and Planetary Science Letters*, 150(3-4), 177-189. doi: 10.1016/
s0012-821x(97)00109-x

Spera, F. (2000). Physical properties of magma. In H. Sigurdsson (Ed.), *Encyclopedia of volcanoes* (pp. 171–190). Academic press.

Stein, A. F., Draxler, R. R., Rolph, G. D., Stunder, B. J., Cohen, M. D., & Ngan, F. (2015). NOAA's hysplit atmospheric transport and dispersion modeling system. *Bulletin of the American Meteorological Society*, *96*(12), 2059–2077. doi: 10.1175/BAMS-D-14-00110.1

Sudhaus, H., & Jónsson, S. (2009). Improved source modelling through combined use of InSAR and GPS under consideration of correlated data errors: Application to the June 2000 Kleifarvatn earthquake, Iceland. *Geophysical Journal International*, *176*(2), 389–404. doi: 10.1111/j.1365-246X.2008.03989.x

Tachikawa, T., Hat, M., Kaku, M., & Iwasaki, A. (2011). Characteristics of ASTER GDEM Version 2. *Geoscience and Remote Sensing Symposium (IGARSS), 2011 IEEE International*, 3657–3660. doi: 10.1109/IGARSS.2011.6050017

Tarantola, A. (1987). *Inverse problem theory: Method for data fitting and model parameter estimation*. Elsevier.

Tobita, M. (2003). Development of SAR interferometry analysis and its application to crustal deformation study. *Journal of the Geodetic Society of Japan*, *49*(1), 1–23. doi: 10.1136/sokuchi1954.49.1

Touloukian, Y., Judd, W., & Roy, R. (1989). *Physical properties of rocks and minearls. data serises on material properties*. (Vols. 1–2). New York: McGraw-Hill.

Tridon, M., Cayol, V., Froger, J. L., Augier, A., & Bachèlery, P. (2016). Inversion of coeval shear and normal stress of Piton de la Fournaise flank displacement. *Journal of Geophysical Research: Solid Earth*, *121*(11), 7846–7866. doi: 10.1002/2016JB013330

Walker, G. P. (1984). Downsag calderas, ring faults, caldera sizes, and incremental caldera growth. *Journal of Geophysical Research*, *89*(B10), 8407–8416. doi: 10.1029/JB089iB10p08407

Wallace, P. J. (2001). Volcanic SO₂ emissions and the abundance and distribution of exsolved gas in magma bodies. *Journal of Volcanology and Geothermal Research*, *108*(1-4), 85–106. doi: 10.1016/S0377-0273(00)00279-1

Walstead, S. T. (1999). *Fractal and wavelet image compression techniques*. Bellingham, Washington: SPIE Optical Engineering Press.

Wauthier, C., Cayol, V., Smets, B., D'Oreye, N., & Kervyn, F. (2015). Magma pathways and their interactions inferred from InSAR and stress modeling at Nyamulagira Volcano, D.R. Congo. *Remote Sensing*, *7*(11), 15179–15202. doi: 10.3390/rs71115179

Wessel, B. (2016). *TanDEM-X Ground Segment - DEM Products Specification Document* (Tech. Rep. No. 3.1). Oberpfaffenhofen, Germany: EOC, DLR. Retrieved from https://tandemx-science.dlr.de/ doi: DOI:10.1002/hyp.3360050103

Wessel, P., Luis, J., Uieda, L., Scharroo, R., Wobbe, F., Smith, W., & Tian, D. (2019). The Generic Mapping Tools version 6. *Geochemistry, Geophysics,*
*Geosystems*, 5556–5564. doi: 10.1029/2019GC008515

Wright, R. (2016). MODVOLC: 14 years of autonomous observations of effusive volcanism from space. *Geological Society, London, Special Publications*, **426**, 23–53. doi: 10.1144/SP426.12

Wright, R., Flynn, L. P., Garbeil, H., Harris, A. J. L., & Pilger, E. (2004). MODVOLC: Near-real-time thermal monitoring of global volcanism. *Journal of Volcanology and Geothermal Research*, **135**(1-2), 29–49. doi: 10.1016/j.jvolgeores.2003.12.008

Wright, T. J., Parsons, B. E., & Lu, Z. (2004). Toward mapping surface deformation in three dimensions using InSAR. *Geophysical Research Letters*, **31**(1), 1–5. doi: 10.1029/2003GL018827

Yun, S., Segall, P., & Zebker, H. (2006). Constraints on magma chamber geometry at Sierra Negra Volcano, Galápagos Islands, based on InSAR observations. *Journal of Volcanology and Geothermal Research*, **150**(1-3), 232–243. doi: 10.1016/j.jvolgeores.2005.07.009

Yun, S. H., Zebker, H., Segall, P., Hooper, A., & Poland, M. (2007). Interferogram formation in the presence of complex and large deformation. *Geophysical Research Letters*, **34**(12), 1–6. doi: 10.1029/2007GL029745

Zebker, H. A., Rosen, P., Goldstein, R. M., Gabriel, A., & Werner, C. (1994). On the derivation of coseismic displacement fields using differential radar interferometry: the Landers earthquake. *Journal of Geophysical Research: Solid Earth*, **99**(B10), 19617–19634.
<figure><img src="image_3.png" /><figcaption>**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.</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>**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. The dike trace represents the upper edge of the quadrangle ("Dtop" in Fig. S7).</figcaption></figure>
a. Inversion 1: Dike

<figure><img src="image_5.png" /></figure>

b. Inversion 2: Dike and fault

<figure><img src="image_6.png" /></figure>

c. Inversion 3: Dike, fault, and reservoir

<figure><img src="image_7.png" /></figure>

**Figure 4.** Opening and shear displacements of the best-fit models for **a.** Inversion 1 (including a single dike), **b.** Inversion 2 (including a dike and fault), and **c.** Inversion 3 (including a dike, fault, and a reservoir). There is no shear displacement on the reservoir, therefore it is not shown in the right-hand column of **c.** The marginal posterior probability density functions can be found in the Supporting Information, Figures S9, S10, and S11.
<figure><img src="image_8.png" /></figure>

**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.
<figure><img src="image_9.png" /><figcaption>**Figure 6.** A simplified schematic illustration summarizing the findings from this study. In 2015, a reservoir (blue spheroid) at depths of 4–5 km b.s.l. fed a dike intrusion and intra-caldera fissure eruption. An extraction of less than 7% of magma from the reservoir is sufficient to un-clamp 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³. 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×.</figcaption></figure>