Research ArticleGEOLOGY

Goldilocks conditions required for earthquakes to trigger basaltic eruptions: Evidence from the 2015 Ambrym eruption

See allHide authors and affiliations

Science Advances  03 Apr 2020:
Vol. 6, no. 14, eaaz5261
DOI: 10.1126/sciadv.aaz5261


Observations indicate a strong correlation between the occurrence of volcanic eruptions and earthquakes. While increased volcanic activity has been observed following both local and distal earthquakes, some of the largest recorded earthquakes aren’t known to have triggered an eruption. Here we investigate whether an eruption and associated dike intrusion at Ambrym volcano was triggered by an Mw 6.4 earthquake which occurred 30 hours earlier. Modeling suggests that stress changes induced by the earthquake were too small to account for the overpressure in the dike without additional bubble growth to pressurize the magma chamber. We find that the magma must be both H2O-saturated and at lower temperatures than those expected for newly intruded basalts. Too hot and the stress drop required to grow the bubbles is too large, too cold and the magma can no longer flow. These observations suggest that partially cooled and crystallized basaltic magmas are more susceptible to triggering from earthquakes.


While the probability of an eruption being triggered by an earthquake is considered low (1), increased numbers of eruptions and unrest episodes following large magnitude earthquakes support short-term triggering at both near-field (tens of kilometers) and far-field (hundreds of kilometers) distances (16). Although the mechanism that ultimately triggers an eruption continues to be debated, perturbations to the stress field are thought to be a primary cause. At long distances, dynamic stress changes due to the passage of seismic waves may induce the nucleation and growth of bubbles within the magma chamber (1, 4) and cause liquefaction of the crystalline mush (7) or the collapse of magmatic bubbles and foams because of the oscillatory motion within the reservoir (8). Although they decay more rapidly with distance, static stress changes are permanent, explaining the time delay in triggering some eruptions (9, 10). Static stressing involves the volumetric expansion of magma chambers (11), promoting the growth of bubbles, the squeezing of magma stored within shallow chambers causing its upward migration and eruption (4), and can also control the location of intrusions (11, 12). While numerous mechanisms have been suggested, for typical stress changes induced by earthquakes, it is thought that the overpressure within the magma body must be close to critical for an eruption to initiate (1). Here, we investigate the possible triggering of a dike intrusion by a local moment magnitude (Mw) 6.4 earthquake, which occurred at Ambrym Volcano, Vanuatu, Southwest Pacific, in 2015.

Ambyrm lies within the Vanuatu volcanic archipelago, an intra-oceanic volcanic arc associated with the eastward subduction of the Indo-Australian plate beneath the Pacific plate (Fig. 1) (13, 14). The island is a 35 km × 50 km wide, predominantly basaltic, shield volcano. A 12-km-wide caldera, located at its summit (Fig. 1), developed ∼2 ka ago (15) during major collapses and explosive events (16, 17) and is now the locus of current activity (18). Over the last century, eruptive activity has been focused around two large cones, Marum and Benbow, located within the caldera (Fig. 1) (16). The cones, which host episodic lava lakes, continuously degas and are cut by multiple past dike and sill intrusions (17). Observations of gas fluxes (18) and seismic data (19) indicate a reservoir located at ∼3- to 4-km depth in the vicinity of Marum and Benbow. The latest eruption of Ambrym occurred in late 2018 when a large dike, ~0.7 km3, was emplaced along the east rift zone and was accompanied by ~3 m of subsidence within the caldera (20). Modeling of the deformation suggested multiple sources from within the caldera fed the intrusion at depths of 2 to 5 km (20).

Fig. 1 Shaded relief map of Ambrym Island and its tectonic setting.

In the main figure, the orange polygons and associated black lines indicate the location of the east and west rift zones. The black line shows the location of the caldera rim, and recently active vents Marum, Benbow, and Niri Taten are shown in red. The yellow dot shows the location of the 19 February 2015 Mw 6.4 earthquake and moment tensor. The map in the top right shows the tectonic setting of Vanuatu and Ambrym (red). The heavy black line shows the location of the plate boundary between the Pacific (PAC) and Australian (AUS) plates. The dashed line shows the location of the Back-arc Thrust Belt (BATB), the light gray polygon shows the location of the incoming D’Entrecasteaux Ridge System (DER), and the arrow shows the plate convergence rate (55).

Before 2018, the last eruption of Ambrym occurred on 20 February 2015 (local date, 21 February). Although it was only classed as a minor eruption, it led to the creation of a new vent within the main caldera floor, which was reactivated in 2018 (20), and caused the first lava flow in over a decade. Around 30 hours before the eruption, an Mw 6.4 earthquake occurred at a shallow depth of ∼10 km south of the volcano. Moment tensors (Fig. 1) suggest a shallow thrust event located on a fault striking ∼N-S, consistent with the geometry of the back-arc thrust belt (Fig. 1).


InSAR data and deformation modeling

Using ascending ALOS-2 synthetic aperture radar (SAR) data acquired on 24 January and 4 March, we form an interferogram (see Materials and Methods and Fig. 2) spanning the eruption on 21 February. Deformation is concentrated along a northwest-southeast trending zone running from the main vents down the southern flank. The asymmetric deformation pattern, consistent with the intrusion of a dike, shows line-of-sight (LOS) displacements of up to 1.5-m south of the main active vents extending along an ~4- to 5-km-long region (Fig. 2). An area of incoherence in the interferometric SAR (InSAR) data along the center of the dike is associated with new lava flows (fig S1). Peak displacements are observed ~2 km south of Marum and 3 km east of Benbow, in the vicinity of the new eruption site. Although the dike rotates toward the rim of the western vent (Benbow), there is no clear evidence of deflation at either vent site.

Fig. 2 ALOS-2 interferogram showing the deformation associated with the dike intrusion on 21 February 2015.

(A) Observed, modeled, and residual. The black lines shows the surface trace of the intrusion, heavy red lines are the location of the main vents (Marum, Benbow, and Niri Taten), and the pink line shows the outline of the lava flow based on the change in radar coherence. (B) Best-fitting distributed opening model for the intrusion. The heavy black lines indicate the locations of the main vent sites. (C) Depth-averaged opening along the length of the dike from north to south.

To model the dike, we first used a Bayesian inversion [Geodetic Bayesian Inversion Software (GBIS) (21); see Materials and Methods and fig. S2] to invert for the position, dip, strike, depth extent, and opening of the dike based on the formulations of Okada (22) (see Materials and Methods). After fixing the geometry, we discretized the dike into ~1 km–by–1 km patches and solved for the best-fitting distributed opening. The model indicates that the dike had a maximum opening of 3.5 m extending from the surface to ~3.5-km depth (Fig. 2 and fig. S2). Less opening is predicted at the northern end of the intrusion where a maximum of 2 m of opening is predicted between 1- and 2.5-km depth. In addition to the dike, a small offset in the phase data extending south from Benbow crater is consistent with shallow slip of 0.5 m on a normal fault. The model explains 96% of the data variance and predicts a total intruded volume of 0.045 km3 and, assuming an average 1-m-thick lava flow with an area of 1.27 km2 based on the amplitude imagery (fig. S1), an additional 0.0013 km3 extruded on the surface.

The region of maximum opening along the dike is offset by ~2 km from Marum and Benbow vents, and neither shows evidence of shallow deflation, suggesting that the magma was fed from below. However, unlike the 2018 eruption, which showed widespread subsidence from across the caldera and at both cones (20), we do not see any indication of deeper deflation feeding the intrusion. To look for the potential source for the dike, we use postdyking deformation derived from ascending Sentinel-1 and descending ALOS-2 data acquired since early 2015 (see Materials and Methods, Fig. 3, and figs. S3 and S4). Two years after the intrusion, both the ascending and descending data indicate a broad region of subsidence of up ~10 cm/year over a ~30-km2 area. Two areas of localized subsidence of up to 20 cm are observed in the first 3 months after the intrusion. The locations are associated with new lava flows near Niri Taten (fig. S4) and an area of 2 km to the southeast (fig. S4) and is likely the result of their cooling and contraction (23). Using the post-eruption InSAR data, we estimated the best-fitting displacement rate (see Materials and Methods) and used the same Bayesian inversion to solve for the best-fitting source based on a horizontal dislocation to represent a contracting sill (Fig. 3). The best-fit model suggests a source at ~3.1-km depth consistent with independent estimates from gas chemistry and seismic observations of the magma storage area (18, 19). The model suggests up to 0.5 m/year of contraction centered in the region of maximum dike opening, south of Benbow (Fig. 3), with lesser amounts of contraction following the trend of the dike intrusion above. The observed post-intrusion subsidence may be related to continued drainage of magma into the dike or cooling and contraction of the residual material left in the sill after feeding the dike. Alternative mechanisms, such as viscoelastic relaxation or magma degassing, could also play a role (24). Regardless of the cause, the coincidence between the location of the inferred sill and base of the modeled dike supports the conclusion that the dike was fed from the same shallow reservoir. Because of the large deformation associated with the shallow dike, it is likely that the shallow intrusion masks the deflation of the deeper magma body. Discrepancies between the volume leaving a magma chamber and entering a dike are commonly observed (2527). Dike volumes can be more than three times the volume predicted from the deflating source (27). In this case, the contraction of a sill at 3-km depth by 0.015 km3 (a third of the volume intruded into the dike) gives a maximum LOS increase of ~20 to 30 cm, only 10 to 15% of the observed LOS decrease caused by the dike. Because of the difference in sign (LOS increase versus decrease), any subsidence from the deeper source will lead to an underestimate in our inferred dike volume.

Fig. 3 Deformation between 2015 and 2017 captured by ascending and descending InSAR data.

(A) Observed ascending displacement rate derived from Sentinel-1 data (see Materials and Methods). (B) Observed descending displacement rate derived from ALOS-2 data. Warm colors indicate motion away from the satellite. The heavy black line shows the outline of the best-fit sill model, and the thin black lines show the location of the lava flow and main vents. (C) Best-fitting distributed sill model explaining the postdiking deformation. Contraction at the edge of the model is a result of the inversion fitting noise outside of the subsiding region. The blue dashed line shows the surface trace of the dike intrusion, and the black lines are the same as in (A) and (B).

Stress interactions

Thirty hours before the eruption, there was an Mw 6.4 earthquake ~15 to 20 km south of the eruption site at ~10-km depth. To investigate the effect of the earthquake on the shallow magmatic system, we calculate the normal stress change along the modeled sill based on the moment tensor estimates from GEOFON (; Fig. 4) and the U.S. Geological Survey (USGS; fig. S5). The earthquake is modeled as a rectangular dislocation in an elastic half space (22) with uniform slip. The fault’s dimensions (13.8 km) are assumed to be equal and are determined using a global magnitude-scaling relationship (28). Assuming a shear modulus of 30 GPa, the slip (~70 cm) is then fixed to conserve magnitude. We test both nodal planes (Fig. 4 and figs. S6 and S7), but based on the tectonic setting (14) and location of the Back arc thrust belt (Fig. 1), the westward dipping nodal plane seems more likely. To account for the uncertainty in hypocenter locations, we systematically move the location and depth of the hypocenter and assess the stress change at for each position (Fig. 4 and fig. S5). Considering a horizontal and vertical location error of ±5 and ± 4 km, respectively, the stress models suggest that there was likely an overall drop in the pressure (increase in normal stress) within the source region rather than a pressure increase with a maximum absolute change in normal stress of between 0.03 and 1 MPa (Fig. 4 and fig. S6). Both the USGS and GEOFON solutions predict similar stress changes (fig. S5), whereas using the eastward dipping nodal plane produces a small reduction in the normal stress of ~0.01 MPa.

Fig. 4 Stress changes associated with the Mw 6.4 earthquake and the effects on bubble growth.

(A) Maximum change in normal stress changes computed across the sill for different hypocenter locations at 11-km depth (others depths and moment tensors are shown in figs. S5 to S7). Each pixel shows the maximum stress change on the sill if the hypocenter was at that at that location. The moment tensor shows the actual location based on the solution from GEOFON and 5-km horizontal uncertainty. (B) Dynamic stress change derived from the USGS ShakeMap peak ground velocity (PGV) values. Contours are at 0.1-MPa intervals. (C) Plot showing the predicted overpressures in the dike based on Eq. 2 for a range of densities and excess magma pressures. The black dashed box shows the range of reasonable densities/density contrasts for a basaltic melt. (D) Crystal volume fraction plotted against temperature for an Ambrym magma composition with different bulk H2O contents. The area between the dashed lines indicates the range of crystal fraction where lock up begins to occur, and the brown box shows the range of temperatures over which the magma remains mobile enough to erupt. (E) Expected volume increase with decreasing pressure for the 2.5 and 1.5 wt % H2O magmas. The colored lines show the effect of different temperature melts. The solid black line shows the inferred total source volume for each pressure drop. Note that for 1.5 wt % H2O, the 950° and 975°C lines plot on top of one another.

In addition to static stress changes, the passage of seismic waves can also generate significant nonpermanent stress changes within volcanic systems (1, 4). Dynamic stress changes decay less rapidly than static stresses (1) and can therefore exceed the static stress change at larger distances. Dynamic stresses can be directly inferred from the peak ground velocities (PGVs) and are equal to (PGV/Vs)μ, where Vs and μ are the shear wave velocity and shear modulus, respectively (29, 30). As there were no PGV observations associated with the earthquake, we use the PGV predicted by the USGS ShakeMap to calculate the magnitude of the dynamic stress change (30). The predicted PGVs in the region of the caldera are 0.09 to 0.15 m/s, and assuming a shear modulus of 30 GPa and a shear wave velocity of ~3000 m/s, the peak dynamic stress change would then be ~0.9 to 1.5 MPa (Fig. 4), consistent with estimates from similarly sized earthquakes at similar distances (30).


For a dike to breakout from a magma body, the excess pressure must be large enough to overcome the tensile strength of the host rock (31). Once the dike has been injected, the maximum opening, b, along an elastic mode I crack as a result of internal fluid overpressure, P0, is given by (32)b=2P0(1ν2)LE(1)where ν, E, and L are the Poisson’s ratio, Young’s modulus, and half-length of the dike, respectively. Assuming a Young’s modulus in the range of 10 to 30 GPa, a Poisson’s ratio of 0.25, a half-length of 2500 m, and a maximum opening of 3.5 m, the final overpressure for the intrusion was ~4 to 11 MPa. To achieve such an overpressure, the excess pressure in a magma chamber is given by (31, 33)Pe=P0(ρrρm)ghσd(2)where Pe is the excess pressure in the magma chamber, ρr is the density of the host-rock, ρm is the density of the magma, g is the acceleration due to gravity, h is the depth extent of the dike, and σd is the differential stress, which is approximately equal to the tensile strength at the surface (34). For a reasonable range of magma densities (2650 to 2750 kg/m3), a rock density of 2800 kg/m3, σd as 1.25 MPa (34), and 3.0 km for the depth extent of the intrusion, the excess pressure in the underlying chamber would need to be in excess of ~2 MPa to achieve the required overpressure (Fig. 4).

For an eruption to initiate, the excess pressure in the magma chamber must be sufficient to rupture the chamber walls. In many cases, this excess pressure is caused by the injection of new material into the magmatic system with deformation and seismic unrest commonly observed before eruptions (35). Alternatively, external triggers such as local or regional earthquakes (16) can generate sufficient overpressures to generate an eruption. Unlike the eruption of Ambrym in late 2018, where there was a large increase in seismicity in the buildup to the eruption (20), there was no precursory unrest before the eruption in 2015 leading the Vanuatu Geohazards Observatory to downgrade the activity level of the volcano. An interferogram from November 2014 and January 2015 also indicates that there was no deformation across the caldera floor (fig. S4). Unfortunately, due to the available SAR acquisitions, we are unable to resolve whether there was any magma buildup in the month immediately preceding the eruption, but the lack of any precursory unrest and the effusive nature of the eruption do not support a rapid influx of material into the system. However, if magma were accumulating at depth, or at low rates over a longer period of time, then it is unlikely that we would be able to detect it with the available data. Although petrological data from past eruptions (18) suggest that magma equilibrates at shallow depths before eruption, it is possible that an intrusion of magma deep in system led to the eruption and earthquake south of the island. The stress change resulting from the deep accumulation of magma has been shown to trigger nearby moderately sized earthquakes (36). To test whether the deep accumulation of magma could have triggered the Mw 6.4 event, which preceded the eruption, we calculate the Coulomb failure stress (CFS) along the inferred fault plane caused by the inflation of a 5 km–by–5 km sill and dike by 0.1 km3 at different locations and depths (fig. S8). In both cases, the general pattern is for the a reduction in the CFS regardless of the emplacement across the region. For the sill, the only instances where we see an overall positive stress change is when the intrusion either intersects or is extremely close to the fault, which places any inflation source well away from the main volcanic center. Similarly, for a dike source, an increase in CFS is only predicted when the intrusion is located adjacent to the fault either at a depth of ~20 km on the hanging wall side or at depths of 2 to 5 km on the footwall side (fig. S8). While we cannot rule out this deeper influx of material from the available observations, based on these analyses, it seems unlikely that an intrusion at depth was responsible for triggering of the nearby earthquake.

If we assume no additional magma intruded deep into the system, then there needs to be an alternate mechanism to generate the excess pressure to drive the eruption. One possibility is the Mw 6.4 earthquake, which occurred 30 hours before the dike intrusion. However, on the basis of our stress analysis, the magnitude of the expected stress changes, either static or dynamic, is too small to generate the inferred excess pressure needed to drive the eruption. Numerical simulations of pressure recovery within a magma body show that the growth of small bubbles following a pressure drop can cause an overpressurization of the magma by an order of magnitude (37). To explore the possible magmatic conditions for this eruption, we use MELTS (38, 39) to model the effect of isothermal decompression of an Ambrym basalt (18) and examine the bulk volume change due to gas exsolution of the magma caused by a decrease in pressure (see Materials and Methods). On the basis of the estimated source depth and melt inclusion pressure estimates (18), we start the decompression models at 100 MPa using a range of initial magma temperatures from 925° to 1025°C with varying water content (see Materials and Methods and Fig. 4). For higher-temperature magmas, (≥1025°C) with 2.5 weight % (wt %) water, gas exsolution does not initiate until the pressure has dropped by 5 MPa (Fig. 4). For the same magma with relatively low water contents, the initial pressure drop required for bubble growth is even larger. At 975°C and 1.5 wt % water, the pressure drop needs to be more than 10 MPa before bubbles would begin to form. However, if the magma is cooler and H2O saturated, then small pressure drops (less than 1 MPa) are sufficient to promote bubble growth, causing an increase in the magma volume and pressurization of the chamber.

As a magma cools, it crystallizes and its ability to flow is reduced (40, 41). Furthermore, the viscosity of magma containing more than ~30 to 40% crystals (particularly tabular feldspars, which are common at Ambrym) will rise significantly to a point of mechanical lock up, which reduces its potential to erupt (42). For magmas with 2 to 2.5 wt % H2O, this places a lower-bound temperature needed for mobility of around 970°C, while for those with closer to 1.5 wt % H2O, the temperature is closer to 1030°C. However, to get bubbles to grow in a magma at this temperature would require a ~40-MPa drop in pressure (see Materials and Methods and tables S1 to S5). Because of the competing requirements for H2O saturation at low temperature and the rheological constraint, our results indicate that there is a very narrow temperature (~50°C) window where magma remains liquid enough to flow and where small pressure drops lead to sufficient bubble formation to pressurize the magma chamber (Fig. 4). Furthermore, H2O-saturated magmas can undergo larger reductions in temperature before locking up, enabling additional bubble growth for small pressure drops (Fig. 4). This may explain why eruptions linked to triggering (3) are typically associated with arc volcanism, whose magmas are generally H2O rich (43), rather than volcanoes associated with ocean ridge systems or oceanic islands whose (mafic) magmas have H2O contents of less than 1% (44). However, note that there are substantially more large earthquakes near arc volcanoes than those in other tectonic settings, making it difficult to assess whether the lack of H2O is a controlling factor.

Another important constraint on whether the intrusion was triggered is the relative timings of the intrusion and earthquake. The surface eruption was reported ~30 hours after the earthquake. Models for the evolution of dike volume indicate that the volume within a dike varies as V(t) = V (1 − e−τ/t) (45), where V is the asymptotic volume and τ is the decay time scale. On the basis of observations of intrusions from Kilauea, Japan, Afar, and Iceland, τ is typically around 4 to 10 hours (45). Assuming the same evolution of volume at Ambrym and decay times of 4 to 7 hours, it would take 20 to 30 hours for the dike to approach its final volume of 0.045 km3 as it intruded toward the surface (Fig. 5). This would allow up to 10 hours for bubbles to form and grow within the magma, consistent with both models of bubble growth in magmas, which suggest time scales of seconds to hours (46), and with the observed time lag between the earthquake and eruption.

Fig. 5 Dike volume evolution with time based on Rivalta et al. (45).

Each colored line is for a different value of τ that controls the evolution time scale. The dashed line shows the time of the eruption after the earthquake, suggesting that after ~30 hours, the dike would have reached its full volume.


Using evidence from the 2015 eruption at Ambrym, we have provided a possible mechanism for triggering of eruptions from small stress perturbations. While our analysis is focused on an Ambrym type magma, these results suggest that fresh magmas intruded into shallow magmatic systems may be too hot to generate significant bubble growth with low pressure drops. Because partially cooled and crystallized magmas will promote bubble growth and subsequent pressurization of a magma body following only a small stress perturbation, they are more susceptible to triggering from earthquakes. In addition, these results indicate that Goldilocks temperature conditions are required, whereby magma remains liquid enough to flow and where small pressure drops can generate enough bubble formation to sufficiently pressurize the magma chamber.


InSAR processing

All of the InSAR data were processed using data using the InSAR Scientific Computing Environment software (47). Topographic corrections were made using the 30 m SRTM version 3 (48). The interferograms were filtered by a power-spectrum filter (49), and unwrapped by SNAPHU (50). To generate the time series (fig. S3), we follow a small baseline approach (51) such that the displacement, at each epoch, is given by(Gκ22)(m0)=(d0)where G is a matrix containing the timespans of each interferogram, ∇2 is a finite difference approximation of the Laplacian smoothing operator, and the factor κ2 determines the weight of smoothing; m is a matrix containing the best fit velocities between each epoch, which, when multiplied by G, gives the displacements contained within the interferograms d. To derive the displacement rate used in the post-eruptive modeling, we used the displacement time series and solve for the best fitting displacement rate at each pixel assuming a constant rate between October 2015 and March 2017 to match the ALOS-2 interferogram. The inversion is weighted using the data variance in a ~2 km–by–2 km region around the reference location outside of the caldera (P4 in fig. S3). For the ALOS-2 data, where we have only a signal long temporal baseline interferogram, we simply divided the interferogram by the total time span.


To model the dike intrusion, we used a Bayesian inversion [GBIS V1.0, (21)] to find the posterior probability distribution of the source parameters. The data were first translated into a local coordinate system and subsampled using a Quadtree algorithm (52). We solved for the position, depth extent, strike, dip, and opening that best fit the observed deformation assuming a Poisson’s ratio of 0.25.

For the distributed opening model, we used the results of the Bayesian inversion to help constrain the geometry of the dike at depth. The surface position of the dike was traced from the offset in the phase data and used the location from the Bayesian inversion as a guide for the central part of the intrusion where coherence is lost. With the geometry fixed, we solved for the best-fitting distributed opening model, m(Axy1κ2000)(mabc)=(d0)where A is a set of matrices representing Green’s functions for the ascending interferogram, which when multiplied by m produce the model displacements at the observation points, x and y. ∇ is the finite difference approximation of the Laplacian operator, which acts to smooth the distribution of slip, the relative importance of which is governed by the size of the scalar smoothing factor κ. a and b are phase gradients in the x- and y-directions, respectively, and c is an offsets to account for the unknown zero phase level. d is a vector containing the observed displacements. We assume a Poisson’s ratio of 0.25 and assign a shear modulus of 30 GPa. For the post-eruption sill, we follow the same method detailed above and discretize the sill into 300 m–by–300 m patches.

Petrological modeling

To constrain both the mineralogy and volatile (predominantly H2O) exsolution of the Ambrym magmatic system, we use Rhyolite-MELTS (38, 39), a thermodynamic model that has been calibrated against equilibrium experimental runs. The vast majority of experiments used for model calibration have been conducted on basalt compositions of which the Ambrym system sits well within (53). The first stage of our petrological modeling requires a representative whole rock composition. For this, we use the composition obtained for a scoria sample, AMB30 from Allard et al. (18), which is a sample of the present activity from Benbow Crater (fig. S6). This is similar to the whole rock composition of an ejected scoria from the 2005 eruption (52). We then use published melt inclusion volatile data to constrain the volatile component of the modeled 2015 magma. Allard et al. (18) report a range in H2O content between ~0.8 and 1.4 wt % from olivine-hosted inclusions, while for the CO2 content of the phenocryst-hosted melt inclusions, they obtain a range between ~500 and 1000 parts per million (ppm). Because of the low solubility of CO2, this range is likely to be an underestimate of the true CO2 content of Ambrym magmas. Nevertheless, CO2 has limited impact on the phase assemblage of the magma, and without further constraints, we chose to use a constant CO2 content of 1000 ppm for all our model runs. Magma oxidation state has a significant role in controlling the final phase assemblage, and for this, we constrain the oxygen fugacity to ΔNNO +1 (1 log unit above the nickel-nickel oxide buffer) based on the olivine-spinel-pyroxene oxthermometer results from Sheehan and Barclay (54). Lastly, the starting pressure for all model runs was at 100 MPa, and because we are investigating the effect of very small pressure decreases, we decompressed at intervals of 1 MPa down to 90 MPa.

Using the above compositional constraints, we then run a series of isothermal decompression sequences at a constant bulk composition and CO2 content while systematically varying the H2O content and the magmatic temperature. It is notable that the modeled liquidus temperature of our magma is ~1200°C using a constrained fO2 of ΔNNO +1. Using this as a starting point, we incrementally decrease the magma temperature in 25°C steps to 925°C. For simplicity, we only show the results of runs from 1025° to 925°C in Fig. 4.

The results of the model are twofold for our purposes. First, we can quantify the exsolved vapor phase and therefore the volume increase because of H2O bubble formation. These data are used to constrain the amount of overpressure exerted by the decompression event. Second, we are able to obtain the mineralogy of the magma at any given condition, including the volume of the phases (crystals and melt). These data are converted into a volume percent, which can then be used to evaluate the physical properties of the magma—the ability for the magma to flow and ascend rather than stagnate.


Supplementary material for this article is available at

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We thank JAXA for access to ALOS-2 data and the EU Copernicus and the European Space Agency for access to and scheduling of Sentinel-1 data. The ALOS-2 original data are copyright of Japan Aerospace Exploration Agency (JAXA) and provided under JAXA RA4 PI Project P1093002. Funding: This work was supported by public research funding from the Government of New Zealand with additional support the Natural Hazards Research Platform to I.J.H. through contract 2015-GNS-02-NHRP and the ECLIPSE Programme, which is funded by the New Zealand Ministry of Business, Innovation and Employment. Author contributions: I.J.H. carried out the InSAR processing and elastic modeling and led the writing of the paper with help from G.K. who modeled the magma composition and rheology. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

Stay Connected to Science Advances

Navigate This Article