From major caldera subsidence to quiescence at the world’s top volcanic degassing source, Ambrym, Vanuatu: the influence of regional tectonics

Eruptive activity shapes volcanic edifices. The formation of 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 multiple intensely-degassing lava lakes within a 12-km wide caldera. Using satellite remote sensing of deformation, gas emissions and thermal anomalies, combined with seismicity and ground observations, we show that in 2018 an intra-caldera eruption preceded 2.5 m of caldera-wide subsidence at Ambrym. Subsidence 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. 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, compressive tectonics can induce rift zone volcanism and subsequent caldera subsidence. Manuscript submitted to Scientific Reports


Introduction
Caldera formation at silicic volcanoes is generally associated with paroxysmal explosive eruptions, in extreme cases producing ignimbrites 1 . For basaltic systems, caldera-forming processes are debated, as geological traces of associated magma discharge are often missing [2][3][4][5] . Ambrym, a remarkably active basaltic volcanic island in the Vanuatu archipelago, hosts two well-defined straight rift zones oriented at N105 • S, radiating bilaterally from a broad 12-km wide caldera [6][7][8][9][10] (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 [11][12][13] .
Continuous open-vent passive degassing from the lava lakes, first reported by Captain James Cook in 1774 14 , ranks Ambrym first in worldwide volcanic sulfur dioxide (SO 2 ) emissions over the past decade [15][16][17][18] . Between 1820 and 1937, 10 extra-caldera rift eruptions were reported at Ambrym 11 . Notably, in 1913, a phreatomagmatic eruption on the island's west coast completely destroyed a large hospital 19 . Based on geochemical evidence, these eruptions are believed to be fed from a central reservoir situated beneath the caldera, and magma transported to the coast by lateral dikes 20 .
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 9 . 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 6 . However, geological evidence for emplacement of voluminous ignimbrites coinciding with this dating is controversial 12 , raising questions about the process of caldera formation. In particular, the relationship between the evolution of Ambrym's caldera, the rift zone's persistence, and the complex tectonic setting of the New Hebrides remains poorly understood.

Results
Precursory intra-caldera eruption On 14 December 2018, a volcano-seismic crisis begins at Ambrym when 8 M<3 seismic events are detected inside the caldera between 13h00 and 20h00 UTC (Fig. 2b). Between 23h20 and 23h40 UTC, Himawari-8 geostationary satellite observations of thermal anomalies and SO 2 emissions indicate the onset of an intra-caldera eruption (Fig. 2c, 2d, S1, S2). Field observations (Fig. 2e) reveal that the eruption initiated along a N110 • pre-existing fracture at 590 meters a.s.l. at Lewolembwi tuff ring (Fig. S3), and was characterized by scoria deposits. Petrological analysis of scoria indicates an erupted magma of basaltic-trachy-andesitic composition (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, 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, 2e, S3).
During this phase, surface deformation is measured with interferometric synthetic aperture radar (InSAR) thanks to the serendipitous acquisition of an image by ALOS-2 satellite at 00h24 UTC on 15 December, about an hour after the eruption's onset 23 (Fig. 2d, Table S1). The deformation field spanning 03 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 (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 (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 17 ). 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, Fig. S3, and an estimated average lava flow thickness of ∼5 m extrapolated from 3D mapping of the previous 2015 intra-caldera lava flow, 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 (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, 3b). Joint inversion of SAR data (Tab. S1) reveals ∼3 m of opening along a >30 km long dike, extending from within the caldera to beyond the eastern coast 23 (Fig. 2a, 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, Fig. 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 south coast of the island, as later confirmed during a field campaign in early February and by a single GPS measurement at Ulei (Fig. 3a, 3c).

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 reaches 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, 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, 3d), 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 24 . 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 depth, with a volume change ranging between 195 and 231×10 6 m 3 (depending on the depth and amount of post-eruptive deformation included in the data, Fig. S7) (Fig. 4c). The 2:1 volume ratio between the dike and deflating Mogi source is in fact consistent with mass conservation assuming standard magma compressibility and host-rock shear modulus 25 . 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, 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 10s of seconds) (Figs. 2b, S14). At least 4 additional LP signals are recorded at station SANVU (∼ 150km away) on 16-17 December. These LP events with CLVD mechanisms are consistent with caldera ring faulting or pressure drop within a reservoir 26 .
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, S7). Exponentially decaying subsidence is consistent with elastic response to magma outflow driven by pressure difference between the central reservoir and the eruption site 27 . 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 28 . 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 (Fig. 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árðarbunga or the 2018 Kilauea eruptions 29,30 . 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.

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 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 ) [31][32][33] . Furthermore, to travel horizontally for long distances, there must be a strong horizontal gradient of driving stress in the direction of magma propagation 34 . 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 ) 35 . When comparing diking events at Icelandic and Hawaiian volcanoes, dikes 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 35 . A similar conclusion was drawn after the dike intrusions and lava lake drainage of the 2002 Nyiragongo flank eruption 36,37 . In spite of enhanced thermal activity and increased lava lake vigor (Fig. S16, S17), an absence of uplift in the months to years prior to the 2018 Ambrym crisis is evidenced (Fig. S9), consistent with this inferred lack of overpressure 38 . 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 17 .
The exceptionally thick 2018 dike 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 39 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). From 14 • 30 S to 17 • 00 S, the DER and WTM perturb the Vanuatu arc, resulting in the late Quaternary uplift of the western Santo-Malekula islands 40,41 . Coevally to this collision, the Pentecost-Maewo island chain was formed by growth of a back-arc thrust belt (BATB) since 1.8 Myr 42-44 . 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 45,46 . 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 (σ 1 oriented N110 • ), revealing the dominant overprint of this background compressive stress 47 . 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 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 48, 49 , caldera subsidence was not associated with a large onland explosive eruption, thereby leaving little geological trace at the surface. Rather, the exceptionally 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 (Fig. S17). Modulation of volcanic activity at Ambrym's lava lakes over historical times 11 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.

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 50 , 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. Fig. 2c is then estimated by calculating the thermal index averaged in two regions corresponding to the lava lakes and lava flows. Fig. S1 shows the same time series of Fig. 2c on a longer time window.

The thermal index time series of
The SO 2 flux proxy is estimated in two steps. First, following 51 , 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: 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 Fig. S2b). The threshold is intended to distinguish the signal associated with the presence of volcanic SO 2 from background oscillations. Fig. S2a shows the SO 2 flux proxy of Fig. 2c 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 Tab. S1. Processing of SAR data from the SLC level to wrapped, unfiltered interferograms is performed using the Interferometric SAR scientific computing environment (ISCE) 52 for ALOS-2 StripMap (SM3) and Cosmo-SkyMed (CSK) StripMap data and Generic Mapping Tool's software GMTSAR 53 for ALOS-2 wideswath (WD1) data, with additional post-processing using the NSBAS software 54 . 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) 55 . Interferograms are filtered with a weighted power spectrum filter 56 , followed by a cascading high-pass filter, especially useful in the areas with a high-gradient fringe rate and on the vegetated flanks 34 . An iterative, coherence-based unwrapping method is then used 34 , 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.

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 57 (See Tab. S2). The forward model includes dislocations 58 and Mogi sources 59 . 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-eruptive 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 (modelling the residual deformation field after removing the synthetic deformation from the deflating Mogi), and the closing of the post-intrusion sill 60 (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 (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 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 x 1 km (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: 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) 61 . 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 postintrusion 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.

GPS and bathymetric data
See Supplementary Methods for more details.

Author contributions statement
TS processed the InSAR data, produced the geodetic models and contributed to writing. RG supervised the geodetic analysis and coordinated writing of the manuscript. MB supervised the thermal and gas data processing and contributed to writing. EG supervised the seismic analysis. YM collected and analyzed petrological data. VB processed the GPS data. FL processed the bathymetric data. MV processed the Geoscope seismic analysis. FD contributed to the InSAR data processing. NH contributed to the Himawari-8 data processing. DT contributed to fieldwork. SC and PL contributed to seismic analysis. BP conducted fieldwork, analyzed coastal uplift and contributed to writing.   Fig. 1c. b. Upper panel: seismicity versus time, as a function of planimetric distance with respect to Marum lava lake. Symbol size increases with earthquake magnitude (dark grey: M>3.5 ; light grey: M<3.5). Green arrow highlights migration of seismicity along the east rift zone. Middle panel: cumulative seismic moment release. Focal mechanisms are from USGS (green) and GCMT (maroon). Lower panel: Absolute value of broadband (green) and low-passed (maroon) seismogram at SANVU Geoscope station (epicentral distance ∼ 150 km). Spikes in the low-passed seismogram indicate detection of long period (LP) events. These events are not included in the cumulative seismic moment curve, which only includes volcano-tectonic (VT) events reported by VMGD. See Fig. S14 for details and station location. c. SO 2 flux proxy (blue) and thermal index of the lava lakes (yellow) and intra-caldera eruption (orange) derived from Himawari-8. d. Satellite images acquired during the course of the eruption. Left: Himawari-8 (multispectral, geostationary). Center: ALOS-2 (ScanSAR). Right: Sentinel-2 (optical). e. Ground observations. Left: gas emissions at Lewolembwi (green star) and lava fountaining associated to vent opening on the SW flank of Marum (blue star). Right: Comparison of lava lake crater at Benbow before and in the final hours of the eruption.   December, a small-volume dike was emplaced in the eastern part of the caldera. The dike breached the surface and fed an eruption lasting ∼24 hours. Deflation of Ambrym's magma plumbing system as a result of the intra-caldera eruption induced a drop of pressure in the magma column, leading to lava lake drainage. On 15-16-17 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. 4c), 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.