From prodigious volcanic degassing to caldera subsidence and quiescence at Ambrym (Vanuatu): the influence of regional tectonics

Eruptive activity shapes volcanic edifices. The formation of broad caldera depressions is often associated with major collapse events, emplacing conspicuous pyroclastic deposits. However, caldera subsidence may also proceed silently by magma withdrawal at depth, more difficult to detect. Ambrym, a basaltic volcanic island, hosts a 12-km wide caldera and several intensely-degassing lava lakes confined to intra-caldera cones. Using satellite remote sensing of deformation, gas emissions and thermal anomalies, combined with seismicity and ground observations, we show that in December 2018 an intra-caldera eruption at Ambrym preceded normal faulting with >2 m of associated uplift along the eastern rift zone and 2.5 m of caldera-wide subsidence. Deformation was caused by lateral migration of >0.4 cubic kilometers of magma into the rift zone, extinguishing the lava lakes, and feeding a submarine eruption in the rift edge. Recurring rifting episodes, favored by stress induced by the D’Entrecasteaux Ridge collision against the New Hebrides arc, lead to progressive subsidence of Ambrym’s caldera and concurrent draining of the lava lakes. Although counterintuitive, convergent margin systems can induce rift zone volcanism and subsequent caldera subsidence.

Broad caldera formation (>10 km in diameter) is often attributed to ignimbrite-forming, explosive eruptions 1 . For mafic to intermediate systems, however, caldera-forming processes may also be linked to the lateral propagation of dikes that arrest at depth [2][3][4][5] . Consequently, geological traces of associated magma discharge are often missing. Understanding caldera-forming processes thus relies heavily on contemporaneous observations, such as in 1968 at Fernandina and 2000 at Miyakejima 6,7 . Ambrym, a remarkably active volcanic island in the Vanuatu archipelago, hosts a 12-km wide caldera, making it one of the largest basaltic shield calderas on Earth 8,9 . Ambrym's caldera has been previously interpreted as resulting from the collapse of a giant tuff cone resulting from a sequence of explosive phreatomagmatic eruptions 10 . Onset of caldera subsidence is dated around 2000 BP, based on two 14 C dates of charcoal embedded in debris flows on the caldera rim and flank 11 . However, geological evidence for emplacement of voluminous ignimbrites coinciding with this dating is controversial [12][13][14] , raising questions about the process of caldera formation.
Ambrym, in addition to its broad caldera, also hosts two well-defined straight rift zones oriented at N105°S, radiating bilaterally 10,11,[15][16][17] (Fig. 1). Ambrym's caldera is the site of intense eruptive activity, with frequent strombolian eruptions originating from at least two permanent lava lakes within or on the flanks of the cones of Marum and Benbow as well as occasional intra-caldera lava flows, most recently in 1986, 1988 and 2015 12,18,19 . Continuous     Fig. S3), and was characterized by scoria deposits associated with lava fountaining. Petrological analysis of scoria indicates an erupted magma of basaltic-trachy-andesitic composition (Supplementary  Table S4). Once the eruption begins, thermal anomalies associated with the lava lakes progressively disappear within 12 hours, suggesting a drop in lava lake level (Fig. 2c, Supplementary Fig. S1). A helicopter flight (on 16 December 03h30 UTC) confirms both drainage of all the lava lakes and the partial collapse of Benbow and Marum (Fig. 2e). A lava flow, accompanied by lava fountaining and producing SO 2 and ash-rich emissions, is also emitted for ~24 hours from a second vent trending nearly N-S at 730 meters a.s.l. on the SE flank of Marum (Fig. 2c,e, Supplementary Fig. S3). During this phase, surface deformation is measured with interferometric synthetic aperture radar (InSAR) thanks to the serendipitous acquisition of an image by the ALOS-2 satellite at 00h24 UTC on 15 December, about an hour after the eruption's onset 27 (Fig. 2d, Supplementary Table S1). The deformation field spanning 3 November to 15 December measures ~1.2 m of motion towards the satellite, consistent with an intra-caldera dike dipping 40°S and a maximum opening of ~2 m, yielding a total volume of intruded magma of ~34 × 10 6 m 3 ( Supplementary Fig. S4). Until then, InSAR measures neither subsidence related to magma reservoir deflation, nor extra-caldera displacement, particularly no motion along the rift zones.
While Sentinel-2 satellite optical images indicate that lava flows reach their full extent on 15 December 23h10 (Fig. 2d), the eruption ends on 16 December around 07h00 UTC with a disappearance of thermal anomalies and an abrupt decrease of SO 2 /ash emissions (Fig. 2c). Spaceborne imaging of SO 2 by Sentinel-5P TROPOMI sensor constrains the total mass of released SO 2 during the eruption to at least ~50-60 kt ( Supplementary Fig. S10), corresponding to a volume of degassed magma of ~13 × 10 6 m 3 (calculated assuming <5% crystal content and 0.075 wt% of sulphur in the melt 23 ). First-order agreement with the volume of erupted material derived from mapping of the flow of ~10 × 10 6 m 3 (constrained by multiplying the spatial extent of ~1.95 × 10 6 m 2 for the lava flow, Supplementary Fig. S3, and an estimated average lava flow thickness of ~5 m extrapolated from 3D mapping of the previous 2015 intra-caldera lava flow, Supplementary Fig. S11) suggests that the magma remaining trapped in the dike did not contribute significantly to the observed degassing.
Triggering of extra-caldera dike injection. Following a M w 5.6 strike-slip earthquake on 15 December at 20h21, a sharp increase in the seismic moment release is detected, marking the beginning of magma propagation into the SE rift zone (Fig. 2b). A few hours before the intra-caldera eruption's end, the lateral propagation of a voluminous dike is evidenced by a migration of seismicity from the caldera toward the eastern tip of the island, reaching 30 km from the eastern caldera rim by 17 December 12h00 UTC. Sub-pixel correlation of a Sentinel-2 optical image acquired on 15 December 23h10 UTC indicates that at that time the dike tip is leaving the caldera but did not reach farther than halfway to the east coast ( Supplementary Fig. S15). On 16-17 December, local inhabitants report progressive fracturing at the coastal village of Pamal (13 km from the caldera border) (Fig. 3a,b). Joint inversion of SAR data (Supplementary Table S1) reveals ~3 m of opening along a >30 km long dike, extending from within the caldera to beyond the eastern coast 27 (Figs. 2a and 4b). These SAR geodetic observations indicate the emplacement of a dike with a total volume of intruded magma between 419 and 532 × 10 6 m 3 (depending on the maximum depth and how far offshore the model extends, Supplementary  Figs. S5, S6). Surface deformation across the trace of the dike is asymmetric, with more deformation to the south (Fig. 3a), indicating that the dike dips ~70° to the south. Due to this asymmetry, coastal uplift in excess of 2 meters occurred along the southeastern coast of the island, as later confirmed during a field campaign in early February and by a single GPS measurement at Ulei (Fig. 3a,c).
Submarine eruption and caldera subsidence. The extremely narrow breadth of the faulted corridor observed above the dike at the surface, as small as 400 m along the east coast (Fig. 3b), indicates that the dike almost reached the surface. However, magma does not erupt from on-shore fractures and only minor gas emissions are detected from space until 17 December 14h00 UTC (Fig. 2c). The end of dike propagation, marked by an abrupt decrease of seismic moment release, takes place around 17 December 16h00 UTC (Fig. 2b). InSAR-derived models predict that maximum opening at the surface occurs offshore (Fig. 4b, Supplementary Fig. S6), suggesting a submarine eruption. This is confirmed on 18-19 December, when basaltic pumice is collected on the beach near Pamal and Ulei (Fig. 3c,d), indicating that lava erupted underwater. Although the depth and exact location of the underwater fissure are uncertain, the nature of erupted material (basaltic pumice) indicates a shallow (<100 m b.s.l) and high-rate underwater magma supply able to sustain a protective gas-rich envelope which allows pumice to cool before contact with sea water, preventing it from sinking 28 . This transport method differs from, for example, floating eruption products collected during basaltic eruptions in the Azores, when gas trapped in hollow cavities caused buoyant basaltic balloons to float to the ocean surface 29 . An alignment of volcanic cones visible in the bathymetry is consistent with an offshore prolongation of the rift zone, suggesting that similar submarine eruptions took place in the past (Fig. 3a).
In addition to the dike intrusion, we also measure >2 m of subsidence at Ambrym's summit craters, consistent with deflation of a nearly symmetrical pressurized source, roughly centered on Marum at ~4.5 km (Japan Meteorological Agency). Center: ALOS-2 (ScanSAR) (JAXA). Right: Sentinel-2 (optical). (e) Ground observations. Left: gas emissions at Lewolembwi (green star) and lava fountaining associated with vent opening on the SE flank of Marum (blue star). Image courtesy of and copyright to John Tasso, Vanuatu Island Experience. Right: Comparison of lava lake crater at Benbow before and in the final hours of the eruption. Maps generated with GMT and edited in Adobe Illustrator. (2019) 9:18868 | https://doi.org/10.1038/s41598-019-55141-7 www.nature.com/scientificreports www.nature.com/scientificreports/  www.nature.com/scientificreports www.nature.com/scientificreports/ magma compressibility and host-rock shear modulus 30 . Furthermore, SAR data indicate that caldera subsidence occurred between 15 December 00h24 and 18 December 06h10, hence was coeval to the extra-caldera dike intrusion (Fig. 4b2). There are several locations to the north of the caldera where the ascending ALOS-2 fringe pattern is discontinuous, indicating ~0.4 m slip along the caldera faults (Fig. 4b1, Supplementary Fig. S8). Two M w ~5.6 earthquakes with vertical compensated-linear-vector-dipole (CLVD) focal mechanisms are recorded in the migration's final 24 hours, on 16 December 15h34 UTC and 17 December 01h49 (Fig. 2b). These events are characterized by a long period (LP) content and long duration (exceeding 10 s of seconds) (Fig. 2b, Supplementary  Fig. S14). At least 4 additional LP signals are recorded at station SANVU (~150 km away) on 16-17 December 31 . These LP events with CLVD mechanisms are consistent with caldera ring faulting or pressure drop within a reservoir 32 .
In the weeks following the dike emplacement, caldera-wide subsidence, reaching ~80 cm at Marum crater and decaying exponentially with a half-life of about 6 days, is measured using ALOS-2 and Sentinel-1 InSAR (Fig. 4c, Supplementary Fig. S7). Exponentially decaying subsidence is consistent with elastic response to magma outflow driven by pressure difference between the central reservoir and the eruption site 33 . Modeling of this subsidence indicates a horizontal sill at ~4.1 km depth deflating with a total volume change of ~85 × 10 6 m 3 (Fig. 4c). This sill-shaped post-intrusion deflation contrasts with the Mogi-shaped co-intrusion deflation, suggesting that the central reservoir consists of several storage levels 34 . Low-magnitude seismicity during this phase is shallow (<6 km depth) and located within the caldera, possibly related to continued caldera faulting visible in post-eruptive interferograms (Supplementary Figs. S12, S13).
Although no additional large-scale deformation is observed along the east rift zone after 22 December, a localized <12 cm discontinuity is measured across the fractures mapped along the SE coast (Fig. 4c), suggesting a continuation of the distal submarine eruption, driving the progressive drainage of the central magma reservoir, similar to the 2014 Bárarbunga or the 2018 Kilauea eruptions 35,36 . Field surveys confirm that the submarine eruption may have continued past the 27 December, as more pumice was observed on 3 February 2019 than on 27 December 2018. At the time of writing (August 2019), there has been no satellite detection of sulphur dioxide above background levels since late January 2019, contrasting starkly with the intense persistent degassing measured over the past decade 24 .

Discussion
The December 2018 Ambrym diking event sheds light on the stress state that prevails at the scale of Ambrym island, while providing insight into the magma storage conditions beneath Ambrym's caldera (Fig. 5). Joint analysis of remote sensing and seismicity demonstrates that the condition initiating the rift intrusion in December 2018 was the creation of an open fracture at Lewolembwi, connecting the caldera's localized magma supply to the rift zone. Once such a connection is made, a blade-like fluid-filled crack (dike) is able to travel tens of kilometers December, a large-volume dike injection took place along the east rift zone, likely fed through the conduit opened by the previous intra-caldera dike. During its propagation and inflation, the voluminous dike intrusion induced intense faulting and fissuring above its upper edge, especially at the coastal village of Pamal, suggesting a shallow depth of emplacement. However, the dike did not produce any onshore eruption nor any substantial degassing. From 18 December, pumice washed up on the eastern shore of Ambrym indicating a submarine eruption. Continued subsidence of Ambrym caldera during the weeks following this eruption, as well as continuation of localized faulting at Pamal (Fig. 3c), suggest a prolongation of the submarine eruption without substantial additional inflation of the main dike, compatible with passive magma transport from Ambrym's central reservoir toward the submarine eruption site. The inset is a zoom on Ambrym's tectonic setting, emphasizing that Ambrym's rift zone orientation is sub-parallel to the regional maximum compressive stress, allowing us to interpret Ambrym as a large tension fracture. The activated plane and focal mechanism of the 1999 thrust earthquake are noted on the sketch in light blue. www.nature.com/scientificreports www.nature.com/scientificreports/ away from its source, neither erupting at the surface, nor degassing significantly to the atmosphere. To sustain lateral magma propagation, pressure in the dike (P m ) must be greater than the host rock's minimum principal stress normal to the dike plane (σ 3 , positive), a difference defined as the driving stress (P d ) (P d = P m − σ 3 ) [37][38][39] . Furthermore, to travel horizontally for long distances, there must be a strong horizontal gradient of driving stress in the direction of magma propagation 40 . Magmatic systems in an extensionally-loaded host rock (low σ 3 ) do not necessarily require high magma overpressures to drive large diking events (high P d ) 41 . When comparing diking events at Icelandic and Hawaiian volcanoes, dike widths tend to be thicker in Iceland (higher P d in Iceland), but dikes reach the surface more often in Hawaii (higher P m in Hawaii). The >3 m dike thickness at Ambrym brings us to the conclusion that P d is high without having a large P m 41 . A similar conclusion was drawn after the dike intrusions and lava lake drainage of the 2002 Nyiragongo flank eruption 42,43 . In spite of enhanced thermal activity and increased lava lake vigor ( Supplementary Figs. S16, S17), an absence of uplift in the months to years prior to the 2018 Ambrym crisis is evidenced (Supplementary Fig. S9), consistent with this inferred lack of overpressure 44 . The fact that the central magma reservoir was able to store a large volume of magma (0.4 km 3 ) and sustain dike propagation over large distances without prior overpressurization is consistent with vigorous supply of a relatively volatile-poor magma 23 .
The thick 2018 dike-similar in width to dikes at plate boundaries in Iceland or Afar 45 -and energetic earthquake migration indicate that lithospheric stresses primed for magma injection drove lateral magma propagation. Assuming a minimum volume of intruded magma in the 2018 eruption (400 × 10 6 m 3 ), we calculate a minimum expansion rate along Ambrym's rift of ~2 cm/yr (assumptions: last event = 81 years ago; SE rift length = 25 km; height = 10 km). Gravitationally-induced extension can be plausibly discarded as a mechanism generating these stresses, due to the lack of curvature of Ambrym's rift zone 46 and to the fact that the 2018 rift zone intrusion led to uplift of the entire south-eastern part of the island, hence working against gravity. On the other hand, tectonic stresses induced by the proximity of regional faults involved in the convergence between the Pacific and Australian plates may provide the stress conditions driving rift development at Ambrym. Ambrym is facing the collision-subduction of the D'Entrecasteaux Ridge (DER) and, farther to the north, of the West Torres Massif (WTM) (Fig. 1a). The DER and WTM perturb the Vanuatu arc, resulting in the late Quaternary uplift of the western Santo-Malekula islands 47,48 and the growth of a back-arc thrust belt (BATB), uplifting the eastern islands of Pentecost-Maewo (Fig. 1a) [49][50][51] . Ambrym's east coast is located at the southern tip of the BATB, which was last activated in the thrusting 1999 M w 7.5 earthquake 52,53 . We note that the focal mechanisms during the 2018 diking event have P-axis orientations consistent with this regional compression (Fig. 2b). This triggered seismicity does not appear consistent with stress change caused by the dike inflation (increased compressive stress oriented perpendicular to the dike, i.e. N20°) but rather reflects the regional stress state (σ 3 oriented N110°), revealing the dominant overprint of this background compressive stress 54 . In this context, Ambrym may be envisioned as a giant tension fracture oriented sub-parallel to the local maximum compression axis σ 1 (inset of Fig. 5). The importance of the regional stress field's maximum compression axis was also emphasised after the 2000 lateral dike intrusion at Miyakejima, likewise situated in a convergent margin setting 55 .
The 2018 diking event illustrates how tectonically-induced stresses drive magma transport into Ambrym's well-defined rift zone, efficiently siphoning magma away from the caldera in a relatively silent manner for an observer at the surface. Similar to past events in hot-spot basaltic volcanoes of the Galapagos, especially in 1968 at Fernandina 56,57 , caldera subsidence was associated with non-explosive activity, thereby leaving little geological trace at the surface. Rather, the elevated rate of Ambrym's volcanic activity witnessed in the past decades contrasts with the near-complete muting of degassing and thermal activities at the surface during and since the December 2018 extra-caldera dike intrusion ( Supplementary Fig. S17). Modulation of volcanic activity at Ambrym's lava lakes over historical times 18 may therefore be reinterpreted as resulting from recurrent pumping of magma into the rift zone, leading to episodic subsidence of the caldera floor. Open-vent degassing in the years prior to the 2018 eruption may have allowed for an increasing availability of non-overpressurized magma, with lava lakes acting as piezometers of a magma plumbing system that is well-connected to the surface. This offers a glimpse into a system where magma ascent, lateral magma transport, and caldera formation are controlled by regional compressive tectonics.

Methods
Himawari processing. Himawari-8 is a meteorological satellite operated by the Japan Meteorological Agency, providing multispectral observations from a geostationary orbit at 140.7°E. The Advanced Himawari Imager (AHI) covers 16 channels spanning visible to thermal infrared, and acquires images every 20 minutes. Images have a resolution of ~2.3 km at the location of Ambrym.
Following 58 , a raw thermal index is calculated by computing the normalized difference between top-of-atmosphere brightness temperatures from infrared channels centered on 10.41 microns and 3.9 microns: In order to mitigate the impact of clouds and diurnal variations in brightness temperature, a background thermal index is estimated by extracting the raw thermal index for a reference pixel situated at the border of the caldera, that is not affected by the thermal anomaly of lava and ash emissions while sharing similar ground properties as the intra-caldera pixels. This background thermal index is subtracted from the raw thermal index, yielding a corrected thermal index. www.nature.com/scientificreports www.nature.com/scientificreports/ The thermal index time series of Fig. 2c is then estimated by calculating the thermal index averaged in two regions corresponding to the lava lakes and lava flows. Supplementary Fig. S1 shows the same time series of Fig. 2c on a longer time window.
The SO 2 flux proxy is estimated in two steps. First, following 59 , a SO 2 column amount proxy is calculated for each pixel by differencing top-of-atmosphere brightness temperatures from infrared channels centered on 10.41 microns and 8.5 microns: .

CA BT BT
(2) SO m m 10 4 8 5 2 Then, the time series of SO 2 flux proxy Q SO 2 in Fig. 2c is calculated by summing, for each acquisition, the value of CA SO 2 for all pixels with a CA SO 2 greater than 4 K within a ~50 × 35 km box centered on Ambrym (dashed polygon in Supplementary Fig. S2b). The threshold is intended to distinguish the signal associated with the presence of volcanic SO 2 from background oscillations. Supplementary Fig. S2a shows the SO 2 flux proxy of Fig. 3c on a longer time window, plotted both on a linear scale and a logarithmic scale.
InSAR processing. The SAR images used in this study are listed in Supplementary Table S1. Processing of SAR data from the SLC level to wrapped, unfiltered interferograms is performed using the Interferometric SAR scientific computing environment (ISCE) 60 for ALOS-2 StripMap (SM3) and Cosmo-SkyMed (CSK) StripMap data and Generic Mapping Tool's software GMTSAR 61 for ALOS-2 wideswath (WD1) data, with additional post-processing using the NSBAS 62 . Topographic fringes are removed using DLR's TanDEM-X 12 meter Global (TDX) DEM (an average of DEMs acquired between 14 January 2011 and 22 November 2014) 63 . Interferograms are filtered with a weighted power spectrum filter 64 , followed by a cascading high-pass filter, especially useful in the areas with a high-gradient fringe rate and on the vegetated flanks 65 . An iterative, coherence-based unwrapping method is then used 65 , which we will call MPD, and is a module in NSBAS (see Supplementary Methods). The final unwrapped interferograms are then geocoded.
CSK interferograms have low coherence across the island, due to vegetation, atmospheric effects, and a high rate of deformation. Therefore, CSK descending pixel offsets, which complement the ascending ALOS-2 measurements, were exploited to measure deformation during the rift zone intrusion and the post-intrusion caldera subsidence. After unwrapping and geocoding, swaths F1 and F2 are merged for ALOS-2 WD1 interferograms. We also merge along-track frames for ALOS-2 SM3 interferograms, to ensure that the far-field signal (Pentecost island to the north and Lopevi, Paama, and Namuka islands to the south) is included in the inversion.
Geodetic modeling. Inversion procedure. We then perform the same inversion procedure for each of the three datasets, respectively corresponding to the (a) intra-caldera dike, (b) rift zone intrusion and caldera subsidence, and (c) post-intrusion caldera subsidence (Fig. 4). To focus on the first-order geometry of the pressure sources, we mask localized, near-field signals that may bias the model misfits (i.e. masking localized deformation within 2 km of the dike trace, deformation close to the craters that may be due to conduit pressurization, or subsidence due to lava flow cooling and compaction).
After downsampling the data using a distance-based averaging, we then run a non-linear inversion to find the first order geometry of the pressure sources 66 (See Supplementary Table S2). The forward model includes dislocations 67 and Mogi sources 68 . We find that the intra-caldera phase (a) is best explained by a single inflating dike, whereas the rift zone intrusion and caldera subsidence in phase (b) require an inflating dike and a deflating Mogi. We find that post-intrusion caldera subsidence in phase (c) cannot be explained by the same deflating Mogi as in phase (b), whereas a laterally extensive sill-like deflating source better reproduces the deformation measured in post-eruptive interferograms.
Following this non-linear inversion, the geometry of the pressurized sources is held fixed, and we then performed a constrained least squares inversion to investigate separately the distributed opening of the intra-caldera dike, extra-caldera rift intrusion (modeling the residual deformation field after removing the synthetic deformation from the deflating Mogi), and the closing of the post-intrusion sill 45,69 (see Supplementary Methods).
During the rift zone intrusion, there was significant deformation offshore, limiting the model's resolution. Therefore, we perturb the geometry of the dike used to model the rift zone intrusion by extending the depth, as well as the length offshore, by several kilometers. We derive four end-member models to determine a reasonable range of volumes, taking into account uncertainties imposed by the lack of model resolution ( Supplementary  Fig. S6). Final volumes range from 419 to 532 × 10 6 m 3 , with an "average" model chosen with a reasonable misfit, which extends 6 km offshore and 6 km along depth.
Once the distributed opening along the rift-zone dike is determined, we perform a final iteration to ensure that the initial removal of the Mogi model's synthetic deformation does not propagate significant model errors into the distributed opening inversion. We subtract the synthetic deformation field derived from the "average" model from the original datasets, and rerun the nonlinear inversion to solve for the Mogi depth and volume change. The volume and depth do not change significantly (<6% and <3%, respectively).
Temporal evolution of post-intrusion subsidence. In addition to an ALOS-2 interferogram and CSK pixel offsets measuring deformation after the rift zone intrusion (Fig. 4c), ascending Track 81 Sentinel-1 images were also acquired every 6 days, starting from 19 December 2018. The magnitude of deformation measured in the 6-day pairs decreases with time, allowing us to investigate the temporal evolution of the subsidence. The full extent of the deformation signal is not captured in some of the interferograms because coherence is limited to within the unvegetated part of the caldera. Conversely, ALOS-2 interferograms exhibit a better coherence, which allows for (2019) 9:18868 | https://doi.org/10.1038/s41598-019-55141-7 www.nature.com/scientificreports www.nature.com/scientificreports/ mapping the deformation outside the caldera. However, the temporal resolution of ALOS-2 data is insufficient to capture the temporal evolution of this transient post-intrusion subsidence signal.
In order to combine the high temporal resolution of Sentinel-1 and the high coherence of ALOS-2, we design a specific strategy by (a) first estimating the parameters (shape, depth) of the deflating source from the deformation pattern visible in ALOS-2 data and (b) solving for the temporal evolution of the deflation source by fitting the deformation signal projected in the line of sight (LOS) against the deformation measured in Sentinel-1 data.
To exploit the high coherence of the ALOS-2 interferogram, we use the same inversion procedure as mentioned above to invert for closing on a horizontal sill, fixed to a depth of 4.1 km, with patches 1 km × 1 km (Supplementary Fig. S4). The depth was derived from an initial non-linear inversion of a closing Okada plane. We then project the synthetic deformation field into the ascending Sentinel-1 LOS, and scale the model to fit the 15 interferograms spanning 19 December 2018 to 18 January 2019. The scaling is determined such that: k k model where k is the interferogram index, u k and u model are the average velocities of the data and model, respectively (i.e. the displacement divided by the time span). Here, γ k is the scalar for interferogram k representing the rate of deformation in the time interval spanned by the interferogram. See Supplementary Methods for more details.
We then multiply each scalar by the total magma volume loss in the initial inversion (−85 × 10 6 m 3 ), and perform a least squares inversion to find volume loss for each acquisition (Fig. 4c) 70 . An exponential is fit to the volume loss vs. time, such that A = 0.144 m 3 , B = 9.55 days, C = 0.1 m 3 , and t is the time is days since 16 December 2018, which corresponds with the inferred time of onset of caldera subsidence based on seismicity. The half-life of the exponential decay is thus ~6.6 days.
Removing post-intrusion subsidence. The datasets spanning the rift zone intrusion and main caldera subsidence also span the beginning of the post-intrusion subsidence. We scale the post-intrusion synthetic deformation field, using the exponential described above, to remove post-intrusion deformation from the co-intrusion ALOS-2 and CSK data. However, after rerunning the non-linear inversion on these corrected interferograms, the RMS was not improved-24.03 (without removal) vs. 24.74 (with removal). This may be because the source geometry does not remain constant throughout the entire post-intrusion subsidence phase (especially in the days immediately after the intrusive event, 18-22 December, which are not covered by the initial post-intrusion ALOS-2/CSK inversion). We therefore decide to proceed without removing the post-intrusion deformation.