High-threat explosive silicic eruptions commonly contain banded pumice, reflecting magma mingling in the conduit prior to or during eruption. Heterogeneities in tuffs have been attributed to the draw-up of compositionally distinct magmas, in which low-viscosity magmas ascend more quickly than high-viscosity magmas. The Rattlesnake Tuff of the High Lava Plains in Oregon (northwestern United States) represents a zoned magma reservoir where at least five different rhyolite compositions are preserved in banded pumice samples in variable mingled combinations. Geochemical gradients recorded across band boundaries in pumice were modeled using a Monte Carlo least-square minimization procedure to find the complementary error function that best fit observed Si and Ba diffusion profiles by iteratively varying the concentration of each plateau (i.e., the concentration on either side of the band boundary), the center and spacing of the diffusion profile, diffusion length scale, and temperature. Modeling indicates maximum time scales between mingling and conduit ascent from minutes to hours. Viscosity calculations for each rhyolite composition confirm that highly viscous rhyolites have longer ascent times than low-viscosity magmas, strongly supporting a model of sequential tapping of a zoned chamber controlled by viscosity.

Explosive silicic eruptions produce pyroclastic density currents that can cover landscapes and threaten human populations and agriculture, and even affect global climate. Despite the risk associated with these events, questions remain regarding the pre-eruptive configuration of magma as well as the timing of magma ascent and eruption. Many explosive deposits contain heterogeneities as (1) chemically or thermally zoned ignimbrites indicating heterogeneities in the pre-eruptive magma configuration itself, and/or (2) the presence of banded pumice indicating the simultaneous tapping of distinct compositions upon eruption (e.g., Lipman, 1967; Hildreth, 1979, 1981; Spera et al., 1986; Grunder and Mahood, 1988; Bacon and Druitt, 1988; Wolff et al., 1990; Hildreth and Fierstein, 2000; Aguirre-Díaz, 2001; Bacon and Lanphere, 2006; Pabst et al., 2008; Huber et al., 2012; Ellis et al., 2014; Bachmann and Huber, 2016). Moreso, the compositional layering (i.e., bands) in banded pumice (which are commonly intricately layered and folded) indicates the fluid mingling of distinct magma types during magma ascent in the conduit itself (e.g., Andrews and Manga, 2014). Furthermore, the preservation of such sharp bands implies that the duration of time for distinct magmas to be juxtaposed was short.

Through mathematical analysis, Blake (1981) attributed zoning in ignimbrites and compositionally banded pumice to the evacuation of magma by a series of successive spherical sampling shells that increase in volume with time and intersect with density-stratified horizontal layers in a magma chamber (Fig. 1A). This way, less-viscous magma draws up more quickly, leading to mingling of distinct magmas in the conduit upon and during ascent and eruption (Blake and Ivey, 1986; Spera et al., 1986). What has yet to be done is to confirm this model (Fig. 1A) by quantitatively determining the time scales of magma mingling and their subsequent ascent, which is important for interpreting petrologic and seismic evidence of magma recharge and eruption triggering (e.g., Eichelberger, 1980; Ruprecht et al., 2008; Wright et al., 2011; Till et al., 2015; Shamloo and Till, 2019; Bardelli et al., 2020; Shamloo et al., 2021).

Diffusion chronometry is a powerful tool for estimating time scales using the chemical gradients recorded in zoned erupted material but is commonly applied to mineral zoning (e.g., Costa et al., 2020, and references therein). Glass zoning (i.e., the contact between distinct bands in pumice), on the other hand, has been underutilized by diffusion chronometry, despite also experiencing elemental diffusion across zone boundaries, which is recorded as chemical gradients rather than an abrupt step in composition. In the case of banded pumice, we assume the chemical gradient represents the timing between the first contact between two rhyolite compositions in the conduit (i.e., magma mingling) and their eruption, when diffusion ceases. Here, we apply diffusion chronometry to the chemical gradients recorded across the boundaries within banded pumice to estimate the timing of magma mingling and ascent. In addition, we test the role of viscosity as proposed by Blake (1981) where less-viscous materials ascend faster than more-viscous materials and lead to heterogeneities in explosive eruptive products. The novel technique of applying diffusion chronometry to vesiculated material poses many additional challenges. Our method serves as a proof of concept, demonstrating the ability to calculate time scales of pre- and syn-eruptive magmatic processes associated with threatening explosive eruptions using a commonly observed rock type.

Spectacular dark-light banding preserved in pumices of the Rattlesnake Tuff make it the target of this study (Fig. 2; Fig. S1 in the Supplemental Material1). The ca. 7 Ma Rattlesnake Tuff eruption is part of basalt-rhyolite volcanism in the High Lava Plains of eastern Oregon (northwestern United States) (Fig. S2; Streck and Grunder, 1995, 1997, 2008; Laib, 2016) resulting in a poorly welded to nonwelded tuff with an estimated magmatic volume of 280 km3. Rhyolites within this single eruptive unit fall into five compositionally distinct rhyolite groups (A–E; Fig. 1B; Fig. S3), characterized mainly by variations in Si, Ba, and Fe but also minor and trace elements including Ti, La, Eu, Ta, Nb, Zr, Hf, Rb, Cs, Th, and U (Fig. S3; Streck and Grunder, 1997; Swenton and Streck, 2022). As many as four groups have been documented in a single sample (Swenton and Streck, 2022). Pumice are finely vesicular throughout (typically ~70% vesicularity), with vesicle sizes ranging from a few microns to >100 µm wide (Fig. S4).

Swenton and Streck (2022) verified that the five distinct rhyolite types (A–E) were likely stratified by density (Streck and Grunder, 1997) and experienced little to no chemical mixing prior to eruption. Their density calculations show that the less-dense rhyolite A, with the most evolved composition, resided at the top of the magma reservoir and rhyolite E at the bottom (Fig. 1A; Streck and Grunder, 1997; Swenton and Streck, 2022). This interpretation is supported by rhyolite A occurring alone at the base of the eruptive unit, indicating it was the first composition to erupt (Streck and Grunder, 1997). The boundaries between bands in rhyolite pumice are sharp to crenulate, suggesting they were not in contact for very long, and mechanical mixing (or mingling) likely occurred within the conduit itself upon and during ascent.

In addition to being density stratified, the pre-eruptive Rattlesnake Tuff magma was thermally zoned (Fig. 1A). Rhyolite E storage temperature was 900 ± 70 °C and rhyolite A was 800 ± 40 °C (rhyolite-MELTS thermometry; Swenton and Streck, 2022). Quartz-albite-orthoclase and clinopyroxene barometry yield pressures between 170 (rhyolites A and B) and 215 MPa (rhyolites C, D, and E), or ~6–8 km depth (Swenton and Streck, 2022). Rhyolite A was likely water saturated, whereas the other rhyolites were barely undersaturated, as inferred from the position of rhyolite compositions relative to the water-saturated minimum in haplogranite ternary space (Streck and Grunder, 1997).

Chemical gradients of major and minor elements were measured across band boundaries using an electron microprobe with variable spacing depending on the available pumice surface (Fig. 2; Table S1 and Fig. S5 in the Supplemental Material). Note that not all measured transects yielded resolvable diffusion profiles but all are included in Table S1. The data in the distal plateaus (i.e., the elemental concentrations on each side of the band boundary) were averaged and assigned a rhyolite composition based on the criteria from previous studies to establish mingled end members and therefore spatial context for the mingling and ascent time calculated (Fig. 1B; Streck and Grunder, 1997; Swenton and Streck, 2022). A full suite of trace elements was not analyzed, and therefore not all the discriminants of previous work are applied.

Diffusion modeling was employed using Si and/or Ba zoning (preferably both where both elemental profiles were resolvable on the probe) across the boundaries between distinct bands within Rattlesnake Tuff pumice samples (Fig. 2). The measured concentration profiles were modeled with a Monte Carlo least-square minimization procedure implemented in Python to find the complementary error function that best fit each observed element profile assuming an initial condition of a step function that physically represents the first contact of mechanical mingling between rhyolites (Fig. 2; Fig. S5). Arrhenius parameters for Ba were used from Magaritz and Hofmann (1978), which were determined for high-silica rhyolites with conditions relevant to the Rattlesnake Tuff (Table S3). Arrhenius parameters for Si were less straightforward, and instead two diffusivities for Si in rhyolite (i.e., Si in rhyolite with 3 versus 6 wt% H2O; Baker, 1991) were tested. Water content for the Rattlesnake Tuff has been placed at 2–4 wt% H2O and no more than 5 wt% H2O based on rhyolite-MELTS modeling (Swenton and Streck, 2022). The exact water content in the conduit between mingled pairs is unknown, but when modeling Si profiles using the diffusivity of Si in rhyolite with 6 wt% H2O, the resulting time scales disagree with Ba time scales. In addition, the general width of the Ba and Si profiles are the same, indicating that they record similar amounts of diffusive relaxation; therefore, we favor the time results using diffusivity of Si in rhyolite with 3 wt% H2O, which yield time scales that agree with Ba time scales.

To best account for sources of error introduced from performing diffusion chronometry on vesiculated material, the model iteratively varied the concentration of each plateau by creating a Monte Carlo synthetic normal distribution of possible plateau values (after Brugman et al., 2022) as well as varying the center of the diffusion profile to account for the uncertainty in distance and the uneven spacing of measurements to avoid vesicles during measurement. A range of temperatures was modeled based on thermometry results for each Rattlesnake Tuff rhyolite group (e.g., 760–970 °C; File S1). The diffusion length scale (i.e., the square root of the product of the diffusion coefficient and time) was varied by calculating a range of diffusion coefficients using the uncertainties in Arrhenius parameters D0 and EA (pre-exponential factor and activation energy). Therefore, rather than reporting a “best-fit time,” the model reported results as a best-fit profile that produced the lowest misfit between measured and modeled gradients from a combination of time–temperature–diffusion coefficient values (Brugman et al., 2022). A cumulative distribution function of times was produced for a given profile, where the 5th and 95th percentile of the distribution was taken as the reported minimum and maximum values in a given time interval (Fig. 3; Table S2).

Lastly, while our modeling approach tries to account for uncertainty in the measured profile itself, we acknowledge the additional limitations to this study (further discussed in File S1). For example, if vesiculation were to have occurred during or after diffusion, then the diffusion length would have been expanded, yielding an apparently longer diffusion time. This is a likely scenario given that the Rattlesnake Tuff magma was stored at relatively shallow pressures and was water saturated. This implies that calculated time scales here represent maxima. Additionally, a sensitivity test was performed by fitting multiple profiles from the same boundary but with variable degrees of vesicularity (File S2). We find time scales within error of each other (e.g., combination CD in Fig. 3), which indicates either fine differences in time scales are lost in the spatial resolution of our measurements, or the differences in time scales are enveloped in the error associated with best-fit time intervals.

In general, diffusion times range from minutes to hours depending on the rhyolite members involved in mingling. The most striking observation is that time scales involving the mingling of rhyolite A (i.e., combinations AE and AD in Fig. 3) are generally longer than time scales involving the mingling of rhyolite B and C (i.e., BE and CD in Fig. 3). For example, the average time scale for mingled rhyolites A and E is graphic min or graphic hr, whereas the average time scale for mingled rhyolites C and D is an order of magnitude shorter at graphic min. Based on its proximity to the conduit, rhyolite A should theoretically have been the first magma to enter the conduit and erupt, which is reflected in the field observations of rhyolite A being at the base of most sections (Streck and Grunder, 1997). Therefore, it is possible that rhyolite A generally spent the most time in the conduit, contributing to the longer time scales observed from our data set.

Another factor that could have led to longer time scales for mingled pairs involving rhyolite A is viscosity—where low-viscosity magmas are drawn up more quickly than high-viscosity magmas as originally proposed by Blake (1981). Viscosity calculations were performed for each distinct rhyolite composition preserved in the banded pumice (Fig. 4; Table S4). As expected, rhyolite A is the most viscous relative to the other rhyolites (i.e., rhyolite A: 105.6–5.8 Pa·s; rhyolite E: 104.6–5.2 Pa·s) based on its chemical composition, water content, inferred temperature, and crystal content (<1 vol%; Streck and Grunder, 1997; Swenton and Streck, 2022) likely contributing to the longer time scales associated with the mingling and ascent of rhyolite A. Additionally, viscosity contrast ratios were calculated for each mingled pair, illustrating that the higher the viscosity contrast, the longer the ascent times. These results potentially provide predictive power to ascent times of magmas based on their viscosities.

Conduit ascent rates are calculated assuming the shallowest pressure recorded by Rattlesnake Tuff pumice marks the roof of the chamber and therefore the base of the conduit (i.e., 170 MPa ≈ 6 km; Swenton and Streck, 2022). Generally, ascent rates range from 101 to 103 m/s (Table S3). These are in broad agreement with ascent rates determined from other explosive eruptions using other diffusion techniques (e.g., Fig. S6; Humphreys et al., 2008; Myers et al., 2016; Newcombe et al., 2020).

We conclude maximum conduit ascent times for an explosive rhyolite eruption are on the order of minutes to hours and are strongly controlled by magma viscosity. While many limitations exist within this approach, this study serves as a proof of concept that diffusion chronometry across compositional boundaries in banded pumice can be used to calculate pre- and syn-eruptive time scales for threatening volcanic eruptions. This work also provides support to numerical models related to magma flow conduits.

1Supplemental Material. Additional figures with discussion of analytical and modeling methods (File S1: Figs. S1–S6), data collected (File S2: Tables S1–S4), and Python script (File S3). Please visit https://doi.org/10.1130/GEOL.S.23749350 to access the supplemental material, and contact editing@geosociety.org with any questions.

The authors thank Vanessa Swenton and Martin Streck for sharing Rattlesnake Tuff samples, Marie Takach for assistance with the electron microprobe, Kevin Pendergast for foundational undergraduate work, Kara Brugman for coding expertise and advice, and Scott Boroughs for collecting additional probe images and data. Thank you to Joseph Boro, Jamshid Moshrefzadeh, and an anonymous reviewer who provided constructive feedback that greatly improved this manuscript. This work was supported by a National Science Foundation (NSF) EAR Postdoctoral Fellowship (1952808) awarded to Shamloo and an NSF grant awarded to Grunder (1939347).

Gold Open Access: This paper is published under the terms of the CC-BY license.