Ascent rates of rhyolitic magma at the onset of three caldera-forming eruptions

The rate at which magma ascends has a strong influence on the manner in which it (eventually) erupts. Slower ascent allows degassing of volatiles from the magma, favoring a more effusive eruption, whereas fast decompression fosters volatile retention and consequently results in more explosive behavior (Eichelberger et al. 1986; Mangan and Sisson 2000; Cashman 2004; Castro and Gardner 2008). Determining the rates at which magma ascends, and how those rates evolve over the course of an eruption, is thus important for understanding eruptive activity and improving monitoring and response for specific volcanoes (Dingwell 1996). Furthermore, the ability to determine ascent rates through the use of erupted materials permits reconstruction of the progression of activity from individual eruptions, including historic events.

Ascent rates have been estimated using experimentally determined rates of breakdown rim formation on hydrous phases, the growth of microlites in matrix melt, and bubble number densities (see review in Rutherford 2008). However, many of these methods are best applied to slower magma ascent rates, hotter systems with lower silica contents, or are heavily influenced by processes in specific regions of the conduit system (e.g., bubbles nucleating around the fragmentation front; Toramaru 2006; Rotella et al. 2014). For large explosive rhyolitic eruptions, ascent timescales are often so short that they remain difficult to constrain with petrological tools. As a result, many studies have used analytical and numerical conduit models to constrain values (5–30 m/s: see reviews by Rutherford 2008; Gonnermann and Manga 2013), or used estimates based on the diffusion rate of H₂O into bubbles (0.7–5 m/s: summarized in Rutherford 2008). Thus, our ability to determine magma ascent rates for explosive rhyolitic eruptions requires the application of a speedometer that can record short timescales and be quenched rapidly after fragmentation.

Several recent studies have exploited volatile exsolution in response to decompression during magma ascent to estimate ascent rates (Liu et al. 2007; Humphreys et al. 2008; Lloyd et al. 2014; Myers et al. 2016; Ferguson et al. 2016). Melt-filled re-entrants (REs; also commonly referred to as embayments) in phenocryst minerals are not sealed by crystal growth and therefore can record late-stage changes in the melt surrounding the host crystal. Such changes of particular interest here are those resulting in gradients in H₂O and CO₂ (and in more mafic examples, S) in the REs that can be modeled to estimate ascent timescales (Liu et al. 2007; Humphreys et al. 2008; Lloyd et al. 2014; Myers et al. 2016; Ferguson et al. 2016). Because pressure-dependent solubilities are well known and precise measurements of H₂O and CO₂ concentrations can be made, modeling of volatile profiles in REs presents a powerful tool for constraining ascent timescales for individual eruptions (Liu et al. 2007).

The use of REs has advanced in the past decade from diffusion modeling of 1–2 H₂O and CO₂ point analyses (Oruanui: Liu et al. 2007; n = 9) and H₂O grayscale calibrated transects using electron microprobe analyses (Mount St. Helens: Humphreys et al. 2008; n = 3), to high-precision nanoSIMS H₂O and CO₂ transects (Fuego; Lloyd et al. 2014, n = 4; Ferguson et al. 2016, n = 4). Here we use concentration profiles for H₂O and CO₂ determined by Fourier transform infrared spectroscopy for 38 re-entrants in quartz and present best-fit diffusion profiles calculated using a decompression model for 31 of them. We focus on samples from the fall deposits of three rhyolitic supereruptions for which extensive fieldwork has been previously conducted, providing a solid framework for integrating calculated ascent rates with inferred eruption dynamics. Our study is motivated by a desire to understand: (1) whether decompression rates increase with increases in eruption intensity inferred from field studies; (2) how ascent rates evolve over the course of an eruption; and (3) whether ascent rates can be related to inferred vent geometry, including the shift into the caldera-forming stages of the eruptions.

We investigated three voluminous, caldera-forming eruptions in this study (Fig. 1):

1. The Bishop Tuff, Long Valley, California (650 km³ magma, 0.767 Ma), where early fall activity graded continuously into climactic eruption (Wilson and Hildreth 1997).

2. The Oruanui eruption, Taupo, New Zealand (530 km³ magma, 25.4 ka), which experienced a time break of months between the first outbreak and subsequent activity (Wilson 2001).

3. The Huckleberry Ridge Tuff, Yellowstone (2500 km³ magma, 2.08 Ma) where the initial activity was prolonged and episodic over days to weeks (Myers et al. 2016; Wilson 2009; Wilson et al. in preparation).

Aside from the contrasts in initial behavior, both the Bishop and Oruanui eruptions exhibit a marked transition, inferred from field evidence, from a single-vent configuration to multiple vents related with caldera formation. Because the corresponding transition in the Huckleberry Ridge Tuff is associated with deposition of hot ignimbrite and has no associated rapidly quenched fall material, no ascent rates could be constrained for the caldera-forming phase of this eruption.

In the Bishop eruption, deposition of fall material and onset into coeval flow activity is inferred to have been a continuous process. The transition from a single vent to multiple vents is documented from lithic evidence (Hildreth and Mahood 1986), with the incoming of fragments from the older Glass Mountain complex near the top of fall unit F8 (Wilson and Hildreth 1997), indicating development of vents around the caldera ring fracture (Fig. 1). However, caldera collapse very likely began earlier, as roughly 2/3 of the erupted volume had already been discharged by this stage (Hildreth and Wilson 2007). Importantly, there is no field evidence for significant depositional breaks (more than a few hours) inferred from any of the Bishop Tuff deposits.

In contrast, the Oruanui supereruption can be subdivided based upon layering in its fall deposits into 10 phases, with time breaks inferred between five of these phases (Wilson 2001). The most significant time break, of weeks to months, lies between the activity of phases one and two, where significant reworking and bioturbation of the earlier fall layer is observed (Wilson 2001). These two phases also involve co-venting of “foreign” biotite-bearing pumices (3–16% of sampled 16–32 mm pumices) sourced from an independent magmatic system 10–15 km away (Allan et al. 2012). The presence of this laterally injected magma suggests that these phases were controlled through rifting processes permitting diking to occur along regional faults (Allan et al. 2012). The transition to caldera formation is piecemeal, however, probably starting in Phase 3, with a marked escalation in the volume, discharge rate, and dispersal power of the eruption.

Deposits of the Huckleberry Ridge eruption are three voluminous ignimbrite members (A, B, and C), with initial pre-A and later pre-C fall deposits (Christiansen 2001). The initial (pre-A) fall deposits studied for this paper preserve evidence in their lower parts for small-scale reworking, suggesting short breaks in deposition during the opening eruption stages (Myers et al. 2016; Wilson 2009; Wilson et al. in preparation). Support for this episodic eruption onset comes also from scatter in measured H₂O concentrations in enclosed melt inclusions, observed in multiple layers through the fall deposit and interpreted to reflect slow ascent (days to weeks) of the first-erupted magma (Myers et al. 2016). In addition, the geochemical clustering of melt inclusion, obsidian pyroclast, and matrix glass compositions throughout the fall deposit suggests that multiple vent systems were simultaneously active and co-depositing tephra, tapping three distinct magma domains at the top of the major magma body (Myers et al. 2016).

Whole pumice clasts were sampled where possible from individual layers (linked to published stratigraphies for the Bishop and Oruanui) through each fall deposit (Fig. 1). Where the grain size was too fine for individual pumices to be sampled (e.g., early Bishop, Huckleberry Ridge), glass-coated loose crystals were separated from samples of bulk fall material from distinct stratigraphic horizons. Nine horizons were sampled in each of the Huckleberry Ridge and Bishop fall deposits, and three in the Oruanui eruption, although not all sampled layers contained usable REs (see details in Supplemental¹ Table 1). Samples were crushed (when necessary), sieved to 500–1000 μm and picked for whole, glass-coated quartz crystals. Quartz crystals were then immersed in isopropyl alcohol to aid in visual inspection under a binocular microscope. Crystals were chosen for this study that contained a RE >100 μm in length that preserved a simple morphology, that is, a rectilinear shape with a wide mouth and lacking internal bends (Fig. 2, Supplemental¹ Fig. 1). Selected crystals were individually mounted in crystal bond and carefully positioned to insure intersection of the entire length of the RE (Figs. 2b and 2c). The crystals were then ground and polished to doubly expose the RE, with only those REs preserving inner and rim glass being used for analysis. Quartz-hosted melt inclusions (MIs) from the same horizons and clasts from which REs were chosen were separately mounted for analysis.

After FTIR analysis, quartz wafers were set in a 1-inch epoxy mount for analysis of major elements using a Cameca SX-100 electron microprobe (EPMA) at the University of Oregon. Operating conditions were 15 kV and 10 nA sample current for Si, Ca, Na, Fe, Al, and K, and 50 nA current for Cl, F, Mg, and Ti. A 10 μm defocused beam was used for all analyses. Na, K, Si, and Al were measured first, and their concentrations were calculated using a time-dependent intensity correction in Probe for Windows (Donovan et al. 2007). Glasses were then analyzed for trace elements by laser ablation-inductively coupled plasma-mass spectrometry (LA-ICP-MS) at Oregon State University using a 50 μm spot size, with four glass standards (GSE-1G, BHVO, ANTO, and BCR) for secondary calibration, ²⁹Si as an internal standard, and GSE-1G as a check standard throughout the run. Re-entrant volatile, major, and trace element data can be found in Supplemental¹ Table 3. The full MI data set, which we use below for comparison with the REs, is presented in Myers (2017).

All quartz-hosted REs from the Bishop and Oruanui deposits (Supplemental¹ Table 3) and Huckleberry Ridge (Myers et al. 2016) are high-silica rhyolite in composition (SiO₂ = 75–77 wt%, volatile-free). For the 38 REs analyzed for H₂O and CO₂, lengths ranged from 100–450 μm, providing 4–17 analyzed spots per RE to define profiles (Fig. 2, Table 1). The RE lengths differed by deposit (Supplemental¹ Fig. 2), with the Bishop (100–400, median 240 μm) and Oruanui (140–450, median 255 μm) containing larger and more variable length REs compared to the Huckleberry Ridge (110–230 μm, median 170 μm). Interior volatile concentrations (H₂O, CO₂) of the REs are as follows: Bishop H₂O = 3.0–5.2 wt%, CO₂ = 60–275 ppm; Oruanui H₂O = 1.9–4.0 wt%, CO₂ = 25–100 ppm; and Huckleberry Ridge H₂O = 2.0–3.4 wt%, CO₂ = 50–450 ppm (Fig. 3). The H₂O and CO₂ profiles from REs (n = 9) in the Huckleberry Ridge initial fall deposits were previously published and discussed (Myers et al. 2016). In this paper, we present refined ascent rates for these REs using the modified decompression code presented below so the results can be directly compared with values inferred for the Bishop and Oruanui samples.

Many REs from Bishop and Oruanui have interior H₂O and CO₂ concentrations that are lower than enclosed MIs in quartz from the same deposits, though there is some overlap (Fig. 3). In contrast, for Huckleberry Ridge, there is complete overlap between RE interior and MI values. However, this is largely a reflection of the unusually wide range of MI H₂O values for Huckleberry Ridge (Myers et al. 2016), in contrast to Bishop and Oruanui, which show much narrower ranges of values for H₂O in MIs. Although enclosed MIs acquire a certain H₂O and CO₂ concentration at the time of trapping, their H₂O concentration can be modified by post-entrapment diffusion of H through the host quartz (Qin et al. 1992; Severs et al. 2007). This can occur during long-term storage and/or slow ascent toward the surface (e.g., Myers et al. 2016). Following Myers et al. (2016), we interpret MIs with lower H₂O concentrations at a given CO₂ content to have been affected by diffusive loss during slow ascent, and we assume that the highest H₂O values reflect the H₂O of the melt inclusions after prolonged magma storage at high temperature. Of the three eruptions studied, the Bishop MIs show the least effects of diffusive H₂O loss during ascent (i.e., the least variation in H₂O for a given CO₂ content), Huckleberry Ridge shows the most, and Oruanui is intermediate (see Fig. 3). Interestingly, only the Bishop fall deposits have interior RE concentrations that extend back to values for the MI H₂O concentrations that likely reflect the storage conditions for MIs before ascent began (Fig. 3).

All measured profiles presented here preserve gradients of H₂O and CO₂ (when present) toward their mouths, with longer REs occasionally showing more variability in their interiors (Supplemental¹ Table 2, Supplemental¹ Fig. 3). There were three REs sampled from the upper portions of the Huckleberry Ridge (sample MM3, n = 1) and Bishop (sample BTMMF9-1, n = 2) fall deposits that have much flatter profiles of H₂O vs. distance, and preserve very low H₂O concentrations (1–1.5 wt%). These REs were sampled from fall deposit layers directly beneath thick ignimbrite (Supplemental¹ Table 1). We interpret these flatter profiles to be the result of post-depositional heating of the fall deposit by ignimbrite (e.g., Wallace et al. 2003), which allowed for continued H₂O diffusion to occur after emplacement (Supplemental¹ Fig. 4). Additionally, there are four profiles (2 from Bishop, 2 from Oruanui) that, although they preserve gradients, have interior concentrations corresponding to pressures below 30 MPa (H₂O concentrations between 1.8–2.1 wt%, no CO₂), suggesting prolonged reequilibration in the conduit at low pressures that approach the depth of fragmentation. The ascent rates calculated from these seven REs are not shown in any of the figures, but their profiles can be found in Supplemental¹ Table 2.

Thirteen of the remaining 31 REs preserve no measureable CO₂. These are all from the early Bishop and Oruanui deposits, in which CO₂ concentrations in co-erupted, fully enclosed MIs are relatively low (maximum CO₂: Bishop F1-F8 120 ppm, Oruanui 200 ppm, Fig. 3) compared to those of the Huckleberry Ridge (maximum CO₂ 800 ppm: Fig. 3). Three REs from the Bishop fall deposit (F9), sampled in a location where no overlying ignimbrite was deposited, have normal H₂O profiles and contain CO₂ contents between 220–280 ppm, higher than all MIs measured from the earlier parts of the fall deposit. These F9 CO₂ contents extend to slightly higher values than are found in F9 melt inclusions but overlap with the high end of the range for coeval Ig2Ea ignimbrite and the low end of the range for late-erupted Bishop Tuff (Wallace et al. 1999; Roberge et al. 2013).

We used a 1D forward diffusion numerical model modified from Liu et al. (2007) for fitting the measured H₂O and CO₂ profiles in each RE. In the Liu et al. model, H₂O and CO₂ diffuse through a melt-filled re-entrant in response to changing external boundary conditions governed by magma decompression. The boundary condition at the contact between the host melt and the rim of the RE is based on the melt H₂O and CO₂ solubility at a given pressure, updated at each decompression step, and assumed to be in equilibrium with the gas bubbles present in the melt outside the crystal (Liu et al. 2007). Equilibrium concentrations are dependent on temperature, pressure, and gas phase composition (Liu et al. 2005). If this assumption were to be relaxed, the dissolved H₂O and CO₂ at the RE rim would likely be more elevated, by an amount that would depend on the degree of vapor-melt disequilibrium. One requirement to ensure that near-equilibrium exchange of volatiles between the rim of the RE and external melt can potentially be maintained is visual confirmation of a bubble wall at the mouth of each RE (Lloyd et al. 2014). Although bubbles were noted at the mouths of Bishop and Oruanui REs (Supplemental¹ Fig. 1), they are less frequently observed in the Huckleberry Ridge REs. However, Myers et al. (2016) showed that using the glass adhering to the quartz rim as the exterior boundary position provided an adequate fit to measured profiles. Diffusion is assumed to be negligible once the sample has reached the fragmentation level, because the timescale between fragmentation and quenching of pyroclasts in the plume is likely to be very short (Liu et al. 2007; Humphreys et al. 2008). Estimates of fragmentation pressure for hydrous rhyolitic magmas are in the range of 10–30 MPa based on the vesicularities of pumice clasts and H₂O contents preserved in matrix glasses (Thomas et al. 1994; Gardner et al. 1996) as well as critical bubble volume fraction (Sparks 1978). Although fragmentation pressure may fluctuate during the course of an eruption (Melnik and Sparks 2002; Dufek and Bergantz 2005), we chose a constant value of 10 MPa and later evaluate how alternative values would affect the inferred ascent rates.

We created a Matlab version (a copy of the code can be requested from M. Myers) of the 1D decompression model of Liu et al. (2007), and updated the diffusivity of CO₂ in rhyolitic melt as a function of temperature, pressure, and water content using Zhang et al. (2007). The updated CO₂ diffusivities result in decompression rates that are up to 1.2× faster than those calculated by Liu et al. (2007) (see Supplemental¹ Fig. 5 and Supplemental¹ Table 4). By comparing our measured H₂O and CO₂ profiles to simulated profiles for various decompression rates and amounts of initial exsolved gas, both considered free parameters in all model runs, we can constrain the ascent conditions that most closely reproduce our FTIR-measured profiles. Because CO₂ diffuses at a slower rate than H₂O, modeling both gradients simultaneously provides more robust constraints on ascent timescales. However, as previously noted, 13 out of the 31 REs measured (all from Bishop and Oruanui deposits) have no measureable CO₂, and were thus modeled based solely on their H₂O profiles. The significance of the absence of CO₂ is discussed further in the discussion section.

Starting pressures determined from volatile solubility relationships are constrained either from co-genetic (i.e., same eruptive layer) melt inclusion H₂O and CO₂ measurements (Myers 2017) or the volatile concentrations in the innermost part of the RE analyzed (see Table 2 and Supplemental¹ Table 2). The RE profiles tend to plateau in the RE interior, suggesting that the plateau values represent the starting H₂O and CO₂ conditions at the time of final ascent. However, to further evaluate the results from both starting conditions (MI vs. RE interior), we model both scenarios. All models were run assuming constant decompression rate and isothermal conditions, with temperature estimates as follows: Bishop 740 °C, Oruanui 780 °C, and Huckleberry Ridge 800 °C (see Supplemental Material¹ for discussion). Although the assumption of isothermal ascent is most probably met by an explosive, plinian eruption (Mastin and Ghiorso 2001), we consider how calculated ascent rates would shift for each system if the magma temperature were 20 °C cooler (to evaluate the effect of cooling on calculated ascent rates).

Using the H₂O and CO₂ concentrations of MIs from each sample as starting conditions requires decompression rates between 0.0002–0.03 MPa/s to reproduce measured RE profiles (Table 2). However, using MIs to define the starting conditions generally yields poorer fits to the data, with only seven RE profiles generating χ² < 1, and four of these yielding comparable fits (within 0.1 χ²) to those produced using the RE interior H₂O and CO₂ as starting conditions. The implications of this result will be discussed later. Given these poor fits, we focus on the best-fit decompression rates that use the RE interiors as starting conditions. Using the H₂O and CO₂ concentrations in the interiors of the REs as starting conditions, best-fit decompression rates for the 31 REs modeled range from 0.003–0.5 MPa/s, with χ² ≤ 1 for all 31 profiles. These decompression rates were converted into ascent rates using the difference between the starting pressure, based on H₂O and CO₂ solubility relationships, and the fragmentation pressure, which was assumed to be 10 MPa. Starting pressures were converted to depth using a mean crustal density of 2600 kg/m³, assuming lithostatic conditions. Using the same conversion, a fragmentation pressure of 10 MPa translates to ∼400 m depth. However, using a lithostatic assumption for converting pressure to depth is likely not valid for estimating the fragmentation depth (Mastin 2005; Gonnermann and Manga 2013). Based on numerical conduit models, the depth of fragmentation ranges between 1 and 2.5 km (Mastin 2005; Gonnermann and Manga 2013; Melnik et al. 2005). Here we adopt the fragmentation depth to be 1 km. It should be noted that the effect of using 1500 m vs. 500 m decreases the distance over which diffusion can occur, which serves to slow calculated ascent rates by a factor of 1.7–2.0×.

Best-fit profiles for all REs yield ascent rates between 0.06–13 m/s, spanning roughly two orders of magnitude (Fig. 5). These values overlap with inferred ascent rates reported in the literature for explosive eruptions across a range of eruptive volumes and magma compositions (Rutherford 2008), but are on the slower end of values based on analytical and numerical modeling for plinian rhyolitic eruptions (5–30 m/s; see reviews by Rutherford 2008; Gonnermann and Manga 2013). The majority (75%) of ascent rates presented here are between 0.4 and 5 m/s. These ascent rates yield magma ascent times as short as tens of minutes to as long as 2.5 h. The fastest ascent rates are found in the upper parts of the Bishop fall deposit, coming from as deep as 7.0 km (160 MPa), whereas the slowest ascent rates (0.06–0.48 m/s) are from the earliest phases of the Oruanui eruption (Fig. 5). Notably, these slowest ascent rates are associated with the initial two phases of the Oruanui eruption, but ascent rates are uniformly higher for phase three of the eruption, indicating that average ascent rates increase upward through the stratigraphy (Fig. 5). The Oruanui phase 3 REs are also strongly correlated with higher RE interior pressures and higher RE rim equilibration pressures (R² = 0.84; Fig. 6). This, however, could weaken with additional measurements (currently n = 3 for Phase 3 REs). For the Bishop and Huckleberry Ridge eruptions there is no strong up-section variation in ascent rate (Fig. 5). For all three eruptions, especially Bishop, the later erupted REs preserve higher interior pressures (Fig. 6; Supplemental¹ Fig. 6).

There is a large overlap in the ascent rates observed in the initial fall deposits from all three systems (Fig. 5). This suggests that although the eruptions had differences in initial behavior inferred from field evidence (start/stop vs. continuous activity), the differences do not correspond to obvious systematic differences in the final ascent rates recorded by REs. However, as noted above, the fastest ascents are recorded in several fall layers from the Bishop Tuff, an eruption that has little evidence for time gaps in deposition. By contrast, the slowest are found in the first two phases of the Oruanui eruption, which had relatively weak eruption styles, were separated by a time break of weeks to months, and have been linked to control by an external rifting event (Wilson 2001; Allan et al. 2012). The third phase of the Oruanui eruption, however, lacks the very slow ascent rates that are observed in the first two phases. The transition from phase two to phase three was associated with a shift from a focused vent area on the northeastern margin of what became the structural caldera to an elongate vent alignment on its eastern edge (Wilson 2001). The association of faster ascents with an inferred escalation in eruption intensity and more extensive vent configuration suggests that final ascent rates may be related to vent geometry and mass flux. However, there is no notable up-section increase in ascent rates observed in the Bishop Tuff, including at the onset of Glass-Mountain-derived obsidian lithics (in upper F8, Fig. 5), which is taken as evidence for caldera ring fracture development. One possibility could be that during this period of the Bishop eruption, fall and flow activity may have been from separate and evolving conduits (Wilson and Hildreth 1997), which would obscure a record of any relationships between ascent rate and changes in conduit and vent configuration. Further investigation of the connection between ascent rate, mass flux, and vent geometry would require application to a system with a simple vent geometry and well-constrained variations in mass discharge rate.

A key result from the RE data is that the melt at the mouths of REs preserves a record of greater apparent quenching depths for REs that experienced higher ascent rates. This is most robust for the Oruanui eruption (R² = 0.84: Fig. 6) and suggests one or both of the following: (1) our assumption of a constant fragmentation pressure (10 MPa) is incorrect, and instead, the fragmentation level became deeper as the eruptions progressed, or (2) for REs that experienced greater ascent rates, there was a greater extent of disequilibrium at the boundary between the external melt and the melt in the mouth of the RE.

The error bars on modeled ascent rates presented in Figures 5, 6, and 7 represent the range of decompression rates that produce acceptable fits to the measured data (χ² < 1), taking into account analytical uncertainties. However, as previously mentioned, there are several factors and assumptions that influence the deduced ascent rates, including: (1) depth of fragmentation (D_f), (2) final pressure of fragmentation (P_f), (3) isothermal ascent (T), (4) starting pressure (P_i), (5) assumption of constant decompression rate, and (6) assumption of vapor-melt equilibrium at the RE mouth. Assumptions 5 and 6 are simplifications we adopted to reduce the number of free parameters. At this stage, we do not attempt to quantify how a more complex model would compare to the simplified model results. Figure 7 summarizes how varying assumptions 1–4 individually for Bishop Tuff REs affects the final ascent rates compared to our preferred model where D_f = 1000 m, P_f = 10 MPa, T = magmatic temperature, and P_i = pressure calculated for re-entrant interior values of dissolved H₂O and CO₂. A decrease in temperature of 20 °C, such as might occur from cooling in the conduit, or a change in the fragmentation depth (500 or 1500 m instead of 1000 m) have relatively small effects on calculated ascent rates, shifting values by a factor of 1.2–2.5×. Using the mouth vapor saturation pressure for the pressure of fragmentation slows ascent rates by a factor of 1.1–3×. Similarly, using the melt inclusions H₂O and CO₂ as starting conditions decreases ascent rates by a factor of 1.1–3.5×. However, it should be noted that our preferred model generally produces better model fits for the measured profiles, especially compared to using MI concentrations as starting conditions. Most importantly, the extent of these variations only serves to shift our reported ascent rates (Fig. 7), but does not reduce the two orders-of-magnitude range that is reflected in the 31 measured H₂O and CO₂ diffusion profiles.

One notable observation from our data set is the significant offset between the starting pressures derived from the innermost H₂O and CO₂ concentrations in the REs and the higher pressures associated with the pre-eruptive storage depth of co-erupted MIs (Figs. 3 and 6). We found that the best fits to the RE volatile concentration profiles were achieved when using the innermost RE concentrations, which tend to plateau to constant values, as the starting condition (Supplemental¹ Fig. 3). In contrast, when the initial H₂O and CO₂ concentrations of MIs were used as starting conditions, the χ² misfits were found to be 1.2–15× worse for all but four REs (Table 2). Although ascent rates that use MI values as starting conditions are consistently slower than those calculated using the RE interior concentrations, and provide poorer fits, all estimates fall within a factor of 6 (Fig. 7; Supplemental¹ Fig. 7).

For the 27 REs where the interior RE concentrations provide better starting conditions for the measured profiles (i.e., produce model profiles that are better fits to the measured data points), a mechanism is then required to explain the shallower starting depths, and associated H₂O and CO₂ concentrations, compared to the values for MIs. We interpret this offset as evidence that there was a deeper, and initially slower, decompression period experienced by the REs, which allowed them to partially or fully re-equilibrate with the external melt prior to final, more rapid ascent (Fig. 8). A similar explanation was called upon to explain the wide H₂O variations measured in MIs within individual fall horizons of the Huckleberry Ridge initial fall deposits (Myers et al. 2016). In Myers et al. (2016), the range of H₂O variations was interpreted to represent variable loss of H by diffusion through the quartz host during slow decompression or shallow storage that led to degassed melt surrounding the quartz. We interpret the offset recorded by values for interior RE plateaus with respect to the co-erupted MIs (the highest values of which represent storage conditions) to result from partial reequilibration as the magma initially fed from the storage region into the conduit system. In this interpretation, the measured RE volatile gradients near the outlets of REs are recorders of the faster, final ascent rate, but do not represent the entire ascent history of a given host crystal (Fig. 8). Note that a two-stage decompression model was similarly required to model the H₂O, CO₂, and S profiles of four mafic REs in olivine from the October 17, 1974, eruption of Volcan de Fuego, Guatemala (Lloyd et al. 2014). Although good fits could be achieved for their H₂O and S profiles using a constant decompression model starting from the MI storage region, to also accurately model the CO₂ profiles, an initial slower decompression period between MI storage and a shallower depth was found to be necessary.

We can estimate the timescale of initial, slower decompression by modeling the time required for the melt in an RE to re-equilibrate from starting H₂O and CO₂ concentrations that are similar to those in the MIs to final values that match the innermost part of each RE (Fig. 8, Table 2). Two approaches were used. In the first, we adopted an instantaneous decompression step followed by a re-equilibration period at constant pressure. In this scenario, the RE starts with uniform H₂O and CO₂ concentrations that are the same as in co-erupted MIs. The exterior boundary condition, following the instantaneous decompression, is fixed to the H₂O and CO₂ concentrations measured in the interior of each individual RE. We consider the RE to have re-equilibrated with the external melt when flat plateaus have developed in the RE interior and where profiles have concentrations of H₂O within 0.1 wt% and CO₂ within 10 ppm of re-equilibration values (i.e., well within the measurement error). This method provides a minimum estimate for the times required for RE re-equilibration, which are <1 to 13 h (Bishop), 4 h to 3 days (Oruanui), and 1.5 to 30 h (Huckleberry Ridge: Figs. 9 and 10).

The second approach is to assume constant slow decompression along a degassing path from the pre-eruptive storage region to the same shallower depth, above which the decompression rate rapidly increases, producing the modeled re-entrant gradient. To model this situation, we applied our constant decompression model to all REs, and used the same starting and ending H₂O and CO₂ conditions and criteria for when reequilibration has been achieved as described above. For Bishop REs, the initial slower decompression rates range from 5.0 × 10^-3 to 4.5 × 10^-4 MPa/s, implying 3 h to 4.5 days of slower ascent (Fig. 10). For Oruanui, the slower decompression rates (1.0 × 10^-3 to 2 × 10^-4 MPa/s) require reequilibration times as long as 1–7 days. Last, the Huckleberry Ridge REs require decompression rates of 5.0 × 10^-3 to 2 × 10^-3 MPa/s, equating to times of 5 h to 1 day, similar to the range preserved using the step function. However, these continuous decompression timescales for the Huckleberry REs could only be calculated for five of nine REs. For the remaining four REs, no plausible degassing path could be found between their “starting MI” state and their preserved RE interior H₂O and CO₂ conditions (Fig. 3). This could be due to significant reorganization of the eruptible melt bodies surrounding the quartz host prior to eruption (see discussion in Myers et al. 2016), meaning that an accurate estimate for their starting melt composition cannot be well constrained. The Bishop and Oruanui REs that require the longest reequilibration time tend to have H₂O and CO₂ concentrations in their interiors that record the shallowest depths and, at least in the Bishop fall deposit, generally lack CO₂ (Figs. 3 and 6). This last observation confirms the requirement for an initially sluggish stage of magma ascent, because CO₂ takes longer to re-equilibrate than H₂O.

For the Bishop Tuff the H₂O and CO₂ concentrations measured in the interiors of REs appear to trend toward greater pressures in samples from higher in the stratigraphy (Figs. 3 and 6). One RE interior concentration from the upper fall layer (F9) even overlaps with the starting MI range (Fig. 6). Interestingly, by this point in the Bishop eruption roughly 2/3 of the total erupted volume had been evacuated. We do not have RE data for such late stages of the Huckleberry Ridge and Oruanui eruptions, as the samples analyzed here represent only the first 1–2% of the total volume erupted. We infer this observation for the Bishop Tuff to represent the transition between where there is the need for a two-stage model to reproduce the RE profiles, to the situation where modeling of the measured profiles can accurately use the MI starting conditions. This changeover is reflected in the diminishing χ² (misfit) from using RE interior volatile concentrations when compared to using the volatile concentrations of co-erupted MIs as starting conditions for F9 REs (Table 2). These changes suggest a maturing of the conduit system, where the initial slow period of magma ascent results from a less developed or less well-integrated conduit system, which evolves over the course of the eruption such that later erupted magma experiences little to no initial slow decompression. This hypothesis is consistent with the qualitative model of Scandone et al. (2007), who argued that there is a development phase, leading up to and during a large explosive eruption, during which the storage region and the conduit system become increasingly interconnected. In the Bishop case, our results show that maturing of the conduit system occurred in parallel with the development of multiple vents.

The scatter in measured H₂O concentrations for fully enclosed MIs in the Huckleberry Ridge Tuff (∼1.0–4.5 wt%) has been interpreted to reflect diffusive losses of H through the quartz lattice during slow ascent and shallow storage on a timescale of days just prior to eruption (Myers et al. 2016). Because diffusive loss of H from MIs occurs on similar timescales as those we have deduced for the initial slow decompression stage experienced by the REs (days), a detailed comparison of the timescales from REs and MIs provides a test of our interpretation regarding initial slow decompression in the three magmatic systems.

Determining the timescale for diffusive loss of H from an MI requires information about the following parameters: (1) estimates for initial MI H₂O values at the time of trapping or following extended storage time at high temperature, (2) estimates of external H₂O concentrations at some lower pressure at which the MI has partially or fully reequilibrated, (3) magmatic temperature, and (4) the size of each inclusion and distance to rim (see Supplemental¹ Material, Myers et al. 2016 and Myers 2017 for details on methods and assumptions used). To determine the initial H₂O concentration for each MI we assumed that H₂O behaved moderately incompatibly through partial loss to a vapor phase during vapor-saturated crystallization (see Myers et al. 2016 for Huckleberry Ridge, and Myers 2017 for Bishop and Oruanui reconstructions). For (2), we assumed that the MIs partially reequilibrated in external melts that had H₂O concentrations similar to the plateau values of H₂O found in the interiors of REs from the same sample. For the Huckleberry Ridge Tuff, 66% (total n = 94) of MIs require residence of >1 day in a partially degassed melt to reproduce their measured H₂O concentrations, with 42 inclusions requiring 1–5 days, and the lowest measured H₂O values requiring up to ∼10 days (Fig. 10). Applying the same model to reproduce the H₂O ranges (∼4.0–6.0 wt%) measured from the first two fall layers of the Bishop Tuff, 65% of inclusions (Total n = 54) experienced no observable diffusive loss, equating to <24 h residence in partially degassed melt in the conduit system (Myers 2017). Of the 19 Bishop MIs that show evidence for diffusive losses, 6 require 1–5 days of lower pressure diffusive loss, and 13 spent >5 days in the conduit system. Last, the first and third fall layers of the Oruanui display a wide range of H₂O contents (∼3.0–5.8 wt%) in MIs for restricted ranges in CO₂, similar to the scatter observed in the Huckleberry Ridge (Fig. 3). Of the 49 inclusions (49%) that experienced >1 day in the conduit, 34 require 1 to 5 days of diffusive loss to a degassed melt, and 15 require >5 days (Myers 2017).

There is considerable overlap in the timescales of diffusive loss recorded by MIs and REs, but for all three eruptions the timescales of diffusive loss for MIs extend to longer times than are calculated for the initial slow decompression stage experiences by the REs (Fig. 10). The general agreement between the different approaches provides strong evidence for an early, slower phase of decompression during the opening stages of these eruptions, which we interpret to result from the lack of fully mature, interconnected pathways from the region of magma storage to the lower parts of the conduit system early in the eruptions. For all three systems, this slow decompression phase occurred over time-scales of hours to days, prior to the final, rapid ascent (<1–3 h) that fed the explosive eruptions (Fig. 10). The extension of the MI values toward longer times than are seen for the REs could have several possible explanations. First, the diffusivity of H in quartz is poorly constrained, and the temperature dependence of the diffusivity is even less well known (see Myers et al. 2016). Second, for Huckleberry Ridge, the initial H₂O and CO₂ conditions for the MIs are less well constrained than for Oruanui and Bishop, because trace element chemistry suggests significant reorganization of multiple magma bodies between the time of MI entrapment and the time of final storage prior to eruption (Myers et al. 2016). Third, REs that experienced a longer slow decompression stage may exist but were not chosen for analysis.

The data set presented here represent results of the first study to determine ascent rates using concentration profiles for H₂O and CO₂ in quartz-hosted re-entrants in rhyolitic magmas, with the aim of constraining the timescales of the early stages of activity for three contrasting explosive eruptions. Our results and comparison with decompression rates estimated from diffusive loss of H from melt inclusions (MIs) demonstrate the potential of modeling gradients in re-entrants to estimate ascent rates and timescales from erupted volcanic products. The majority of reentrants (75% of n = 31) require ascent rates between 0.4 and 5 m/s, equating to ascent times in the conduit of tens of minutes up to a few hours. Re-equilibration of the interiors of re-entrants (REs) to lower H₂O and CO₂ conditions when compared to their co-erupted MIs suggests, however, that the majority of REs experienced an initially slower ascent period, on the order of several hours to several days, prior to final ascent and quenching. This inference is supported by the scatter in measured H₂O concentrations from enclosed MIs, also interpreted to reflect slow initial decompression accompanying ascent from the magma reservoir (Myers et al. 2016; Myers 2017). Ascent times estimated solely by using the measured volatile gradients preserved in the REs are thus minimum values in these systems, and can miss hours to days of slower initial ascent conditions.

By collecting samples from stratigraphically controlled levels in the deposits we are able to compare our modeled ascent rates with field interpretations of the progress of each eruption, to provide context for our results. We observe that in the Oruanui samples there is a noted increase in the ascent rates for the layer deposited during periods of increased eruptive flux (Phase 3). The faster ascent rates are also correlated with greater inferred pressures in the interiors and mouths of the REs. The slowest ascent rates determined in this study come from the first two Oruanui eruption phases, which were interpreted to have been mobilized by an external rifting event (Wilson 2001; Allan et al. 2012), and as such, may have involved magma that was not strongly overpressured. The fastest rates modeled in this study (8–13 m/s) are from the Bishop Tuff, the only system for which our sample set covers much of the eruption duration and the only of the three eruptions that shows no field evidence for long time gaps in deposition (Fig. 5). Additionally, late-erupted Bishop REs (fall layers 8 and 9, close to the time of caldera formation) have interior H₂O and CO₂ concentrations that approach pre-eruption magmatic values. We interpret the temporal decrease in the timescale of the initial slow decompression in the Bishop eruption to represent the maturing and increasing interconnectivity of the conduit system. Ascent rates derived from modeling of re-entrants thus provide key insights into the behavior of magma in the hours to days before it reaches the surface. Linkages between RE and MI timescales and field evidence offer promise in understanding the dynamic processes involved in the initiation of voluminous eruptions, as well as providing first-order estimates of the timescales of unrest that might be used for improving warning of impending activity.

We thank Christy Hendrix and Stacey Gunther in the Yellowstone Research Office for permission (YELL-05248) to work in Yellowstone National Park. Financial support was provided by National Science Foundation grant EAR-1524824 to P.W. The presentation of this manuscript was greatly improved with the constructive and thorough reviews by J. Hammer, H. Gonnermann, and an anonymous reviewer.