Explosive-effusive-explosive: The role of magma ascent rates and paths in modulating caldera eruptions

Magma storage depth, pre-eruptive volatile budgets, and magma ascent rates are known to be controls on eruption explosivity (Burgisser and Degruyter, 2015). These three parameters can be interdependent; for example, the abundance of volatiles dissolved in the melt determines the pressure and extent of volatile exsolution during magma ascent, which influence magma buoyancy and ascent rate (Gonnermann and Manga, 2012). Estimating the three parameters is challenging because none of them is measurable and they are typically calculated by combining geophysical, petrological, and/or geochemical methods (e.g., Endo and Murray, 1991; Rutherford, 2008; Armienti et al., 2013) that mostly require lengthy analyses and several corrections and/or assumptions (Edmonds et al., 2019). Thus, it is desirable to find additional tools that can constrain all the above-mentioned parameters.

Magma volatile contents and storage depths are traditionally studied via melt-inclusion analysis. Recently the volatile-bearing mineral apatite has been found to be a good proxy not only for melt volatile budgets (Scott et al., 2015; Stock et al., 2018; Li and Costa, 2020) but also for magma ascent rates (when combined with zoning-pattern analysis; Li et al., 2020a, 2020b). We investigated the volatile chemistry of apatite from the 2006 CE eruption at Rabaul (Papua New Guinea; 4.2004°S, 152.163°E). This eruption was one of the most explosive events in a 20-year-long episode starting from 1994 (Global Volcanism Program, 2006) and had three phases: (1) a sub-Plinian onset that lasted for 12 hr, (2) a ~12 h mixed Strombolian-effusive phase, and (3) discrete Vulcanian explosions over two weeks. The eruption sequence was well documented by scientists at the Rabaul Volcano Observatory (RVO) and has been studied to reconstruct the volcano plumbing system and pre-eruptive magmatic processes (Bouvet de Maisonneuve et al., 2015). The presence of a large dacitic reservoir located at 3–6 km depth with magmas at ~950–1000 °C, which is frequently recharged with basaltic magmas (Itikarai, 2008; Johnson et al., 2010; Patia et al., 2017), was proposed. The minor compositional differences between products from the three phases led Bernard and Bouvet de Maisonneuve (2020) to propose that slight differences in the initial phenocryst abundances and melt water contents in the reservoir are responsible for the change in eruption styles. We investigated the cause of the eruption dynamics by quantifying magma storage conditions, time scales, and ascent rates using apatite.

We studied samples that are representative of the sub-Plinian, effusive, and Vulcanian phases of the 2006 eruption: pumices in a tephra fall deposit, fragments of a lava flow, and fragments of a ballistic cowpie bomb (see a description in Bernard and Bouvet de Maisonneuve, 2020). Major element and volatile (including Cl and F) concentrations in matrix- and clinopyroxene (Cpx)–hosted apatite crystals, melt inclusions, and matrix glass (this study; Bouvet de Maisonneuve et al., 2015) were determined by an electron probe microanalyzer (EPMA). In total, 218 single-point analyses (with 10 μm beam size) were acquired from 161 apatite crystals (40–64 crystals per rock sample), which appear as either inclusions in Cpx or matrix-hosted microphenocrysts (>30 μm in length). Concentrations of OH were calculated from apatite stoichiometry (see the Supplemental Material¹). These results were used to calculate melt H₂O contents using the ApThermo model of Li and Costa (2020). Wavelength-dispersive spectrometry (WDS) X-ray maps of F and Cl (Figs. 1A–1D) and cathodoluminescence images were used to select zoned crystals for acquiring concentration profiles (1 μm beam size, nearly parallel to the crystal c-axes). We used the one-dimensional multi-component (F-Cl-OH) diffusion model ApTimer (Li et al., 2020a) to fit the measured profiles (n = 16) and to estimate magma degassing time and ascent rates. Crystallographic orientations of apatite were accounted for in diffusion modeling (due to strong anisotropy of Cl diffusion; Li et al., 2020a), and root-mean-square deviation (RMSD) was used to assess the goodness of fit. Additional details on the samples and analytical methods are provided in the Supplemental Material.

Apatite crystals in Cpx from the pumice (sub-Plinian), bomb (Vulcanian), and lava flow (effusive) show a similar and narrow range of F-Cl-OH contents (Figs. 1E and 1F) but differ in zoning patterns. Matrix-hosted apatite in the lava shows decreasing OH and increasing F toward the rim of individual crystals (Fig. 1F). WDS X-ray maps show a decrease in Cl concentration from crystal core to rim, followed by a sharp and narrow increase at the rim, most visible parallel to the c-axis (Fig. 1B). The Cpx- and matrix-hosted apatite in the pumice and bomb are not zoned in F or Cl, based on X-ray maps (Figs. 1A and 1C). Concentration profiles in the bomb's matrix apatite show narrow (~3 μm wide) zonation at crystal rims, which was not observed in any crystals of the pumice.

The pre-eruptive melt H₂O contents calculated from apatite (Fig. 2) range from 1.3 (± 0.5) to 3.0 (± 1.1) wt%, except those from the lava's matrix apatite rims at 0.4 (± 0.2) wt%. These agree with measurements in plagioclase- and pyroxene-hosted melt inclusions and in matrix glass in the pumice (Bouvet de Maisonneuve et al., 2015). Except for the lava, the Cpx- and matrix-hosted apatite of the sub-Plinian and Vulcanian eruptions give similar melt H₂O contents, implying limited extent of water exsolution before the two events, in contrast to the effusive one.

With evidence of roughly 10 × longer zoning distances of F-Cl than rare earth elements in our crystals (see details in the Supplemental Material), we considered diffusion (rather than growth) as the main contributor to the observed F-Cl zoning. Rimward F-Cl-OH profiles of the lava and bomb matrix apatite were fit by diffusion modeling to estimate time scales that can be linked to magma degassing during ascent (Fig. 3). The Cl profiles of matrix apatite (Fig. 3) show a “trough” ~3–10 μm away from the crystal rim, which can be explained by the complex coupled diffusion of F-Cl-OH (Li et al., 2020a). Using constant boundary conditions in our model, we obtained well-fitted Cl profiles for apatite from the bomb but poorer fits for the lava apatite. This potentially implies varying boundary concentrations with time during magma degassing and requires future modeling efforts. Using 975 °C (two-pyroxene thermometry), we obtain diffusion times of ~9–28 h for the bomb samples (n = 6) and much longer times for the lava samples (~650–2200 h ≈ 1–3 months; n = 6). The absence of zoning in the pumice apatite implies much shorter times: at maximum 8 h, assuming the zoning distance measurable from this study is 2 μm (n = 4). The main source of error in the time scales comes from the uncertainty in temperature (T; Fig. 4): decreasing T by 25 °C increases the time estimate by a factor of 1.5–2 (e.g., ~3–10 months at 925 °C for the lava).

Water solubility in silicate melts is mostly controlled by pressure and melt composition (Newman and Lowenstern, 2002); the latter plays a negligible role in this study because all deposits have similar bulk-rock and glass compositions (Bernard and Bouvet de Maisonneuve, 2020). Using the average H₂O contents calculated from Cpx-hosted apatite core compositions for the three eruptive phases, we can estimate volatile saturation pressures using the H₂O-CO₂ solubility model of Papale et al. (2006). At 975 °C and assuming no CO₂ in the melt, we obtained H₂O saturation pressures of 30–60 MPa for all three deposits. The addition of ~150 ppm CO₂ in the melt (the mean of CO₂ in pumice Cpx-hosted melt inclusions; Bouvet de Maisonneuve et al., 2015) increases volatile saturation pressures to 50–100 MPa, corresponding to 2–4 km below the crater (taking a mean crustal density of 2600 kg/m³). These agree with the 3–6 km depth of the main reservoir at Rabaul derived from geophysical and other petrological tools (Saunders, 2001; Itikarai, 2008; Bouvet de Maisonneuve et al., 2015). Melt H₂O contents derived from the lava's matrix apatite rims (0.4 ± 0.2 to 1.4 ± 0.5 wt% H₂O; Fig. 2) give us H₂O saturation pressures that translate to magma storage depths of 0.2–1.5 km. These are much shallower than those of the more explosive phases but are still beneath the surface. This implies that the magma that fed the lava flow resided in a storage zone close to the surface where the melt degassed and that such storage did not happen before the more explosive phases.

The magma storage depth and time scales of degassing and residence we reported above provide insights for better interpreting the monitoring data of Rabaul caldera, which is continuously monitored by RVO (McKee et al., 2018). Prior to the 2006 eruption, an uplift of ~15 cm was recorded in January–September 2005 and February–October 2006 (GPS time series in Fig. 4). Such deformation was likely related to mafic magma replenishment up to 10 days prior to the eruption (Fig. 4; Bouvet de Maisonneuve et al. 2015). On a shorter time scale, a threefold increase in the real-time seismic amplitude measurement (RSAM; Fig. 4) started 8 h before the sub-Plinian onset of the eruption and likely reflects fast magma movement to the surface. The time scale is in accord with our time estimates (<8 h) from the unzoned pumice apatite (Fig. 4B).

Following the initial explosive sub-Plinian phase, a lava flow was erupted from the same vent (similarly to Cordón Caulle, Chile; Castro et al., 2013) in a span of a few hours (Global Volcanism Program, 2006). From matrix-hosted apatite, we found evidence of melt water degassing for 1–10 months prior to eruption at 0.2–1.5 km depth (Fig. 4). Given the deeper original depth found for this magma (2–4 km, according to Cpx-hosted apatite), this timing implies that the effusively erupted magma started ascending much earlier than the sub-Plinian magma batch. The spatial resolution of our estimates does not allow discrimination between a slow, multi-step ascent versus a single-step ascent to intermediate depths. Nevertheless, magma movement and storage at these depths for several months allowed further degassing and crystallization, which eventually resulted in an effusive eruption. The longest degassing time of ~10 months suggests magma ascent could have started by the end of the first uplift period between January and November 2005, when Vulcanian activity and ash venting occurred (Fig. 4A). The subsequent rise of the volatile-rich and buoyant sub-Plinian magma batch bypassed this more degassed magma and may have opened and/or widened the pathway, allowing the effusive magma batch to ascend more easily (Fig. 4C; e.g., Thomas and Neuberg, 2014).

Following the effusive phase, Vulcanian eruptions occurred with decreasing frequency over two weeks (Fig. 4D). The melt H₂O contents calculated from Cpx-hosted apatite in the bomb indicate similar depths of pre-eruptive storage as those of magmas of the other eruption phases. Instead of stalling at 0.2–1.5 km depth, this magma was able to ascend to the ground surface within a day, likely because the system had become more open after the ascent of the sub-Plinian and effusive magmas.

Using 3–6 km as the reservoir depth and ≤8 h of ascent time for the sub-Plinian magma batch, we estimate a minimum ascent rate of 0.1–0.2 m/s. Bernard and Bouvet de Maisonneuve (2020) found an ascent rate two orders of magnitude higher of ~40 m/s using permeability and clast sizes (Rust and Cashman, 2011). The two estimates are not inconsistent because the latter one was calculated at the fragmentation depth where ascent velocity is the greatest whereas the diffusion model yields a maximum time (given no zoning in the sub-Plinian apatite) and thus a minimum ascent rate. For the Vulcanian phase, the 9–28 h of magma ascent translate into an average velocity of 0.03–0.1 m/s, in agreement with that of 0.1 m/s found by Bernard and Bouvet de Maisonneuve (2020) from microlites. The drop of ascent rate from the sub-Plinian to the Vulcanian magma batch, despite their similar pre-eruptive melt volatile contents, could be explained by a drop in overpressure in the system after the sub-Plinian event (intense deflation shown by GPS data; Fig. 4), with possibly significant changes in the system configuration such as partial collapse and/or closure of the conduit (de' Michieli Vitturi et al., 2010; Aravena et al., 2018). For the effusive phase, the long zoning distances in apatite include those of the final, decompression-related zoning caused by magma degassing from intermediate depths to the surface. The difference in ascent rates of the eruption phases, resulting in drastic changes in eruption dynamics, can be attributed to a combination of (1) different depths of pre-eruptive residence, which allowed parental magmas to lose volatiles at different rates (sub-Plinian and Vulcanian versus effusive), and (2) a drop in mass eruption rate due to decrease in overpressure (sub-Plinian versus Vulcanian). This is fundamental for interpreting monitoring data at calderas in terms of pre-eruptive magma transfer, especially in the case of extended periods of unrest.

We have shown that the eruption sequence of Rabaul caldera at the surface, from explosive (sub-Plinian) to effusive (lava) and back to explosive (Vulcanian), is distinct from, and thus does not reflect, the sequence of magma movement at depth. We show evidence of an early ascent and stalling at intermediate depth of the magma that led to the effusive lava eruption; this was followed several months later by the ascent of the sub-Plinian and Vulcanian magma batches. Our finding has important implications for interpretation of monitoring data in terms of subsurface magma-transfer dynamics. Lava effusion following the sub-Plinian eruption was the result of several months of magma storage and degassing, the timing of which is consistent with deformation changes, potentially suggesting a connection to mafic magma recharging the system during that period. The fast ascent of the sub-Plinian magma (within a few hours) agrees with the abrupt escalation of real-time seismic amplitude measurement data. Our study also shows that using commonly accessible instruments (scanning electron microscopy and EPMA) combined with thermodynamic and diffusion modeling of apatite, we can obtain volatile contents and degassing time scales of magmas, which are very difficult to acquire simultaneously otherwise.

We thank M. Streck and M. Cassidy for insightful reviews, and M. Norman for editorial handling. This work was supported by the National Research Foundation (NRF) Singapore and the Singapore Ministry of Education under the Research Centres of Excellence initiative, as well as the NRF grant NRF-NRFF2016-04. This is EOS contribution number 440.