Pre-eruptive rhyolite magma ascent rate is rapid and independent of eruption size: a case study from Ōkataina Volcanic Centre, Aotearoa New Zealand

Volatile measurements in mineral-hosted sealed melt inclusions, and open-ended embayments, have previously been used to study magma ascent dynamics in large rhyolitic eruptions. However, despite occurring more frequently, smaller-volume explosive events remain under-studied. We present magmatic volatile data from quartz-hosted melt inclusions and embayments for eight post-25.4 ka rhyolitic eruptions at Ōkataina Volcanic Centre, Aotearoa New Zealand. Seven originated from within the main caldera, and the other erupted from the associated Ōkareka Structural Embayment. Melt inclusions preserve volatile contents of 2.92–5.82 wt% H2O and 13–126 ppm CO2, indicating pre-eruptive storage depths of 4.5–7.4 km, with younger eruptions being more shallow. The lack of correlation between H2O, CO2, inclusion size or distance to the crystal rim suggests magma bodies experienced variable degrees of degassing during magma storage, with some amount of post-entrapment volatile modification prior to and concurrent with final magma ascent. Diffusion modelling of measured H2O gradients in melt embayments indicates ascent rates of 0.10–1.67 m.s−1 over time spans of 20–230 min for the intra-caldera events. In contrast, ascent rates for the eruption from the Ōkareka Structural Embayment may be more rapid, at 1.59–4.4 m.s−1 over a time span of 22–34 min. Our findings imply that the final, pre-eruptive magma movement towards the surface could be less than a few hours. Comparisons with published data for caldera-forming explosive events reveal no clear relationships between final ascent rate, eruption size or initial volatile content, implying that other factors besides eruption volume control rhyolite magma ascent.


Introduction
Studies of volatile concentrations (H 2 O, CO 2 , Cl, F and S) in melt inclusions and volcanic glasses provide useful information about many aspects of magmatic systems (e.g. Gonnerman and Manga 2013;Edmonds and Wallace 2017). Such information includes constraints on magma storage depths and patterns of magma ascent (e.g. Myers et al. 2018Myers et al. , 2019Saalfeld et al. 2022), degassing trends, reservoir dynamics and eruptive behaviour (Blake 1984;Hilton et al. 2000;Lowenstern 2000;Hartley et al. 2014;Chiodini et al. 2016;Grocke et al. 2017;Mittal and Richards 2019;Geshi et al. 2021;Lerner et al. 2021;Waelkens et al. 2022) and volatile emissions prior to and during past eruptions (e.g. Barclay et al. 1996;Saito et al. 2001;Poland et al. 2009;Johnson et al. 2011;Bégué et al. 2014). Determination of magmatic volatile contents also provides insights into the likely hazards associated with volcanic activity, as volatiles are one of the key controls on eruption style and explosivity (Sparks 1978;Blake 1984;Dingwell 1996;Cassidy et al. 2018;Popa et al. 2020Popa et al. , 2021. A key component in linking magmatic volatile abundances to eruption style lies in the use of disequilibrium Editorial responsibility: J.E. Hammer features to measure magma decompression and ascent rates at various stages between the upper-crustal magma storage depths and the surface. Most methods used for assessing magma decompression rates can only record either slow decompression conditions, for example disequilibrium textures such as breakdown rims in mineral grains or microlite growth from the melt (e.g. Rutherford and Hill 1993;Rutherford and Devine 2003;Castro and Gardner 2008), or the very rapid, immediately pre-fragmentation portions of the magma ascent, such as erupted pyroclast textures or bubble number density experiments (e.g. Toramaru 2006;Giachetti et al. 2010). These methods are unable to capture the more intermediate rates of magma ascent occurring during the main part of the magma rise journey through the conduit prior to the more shallow region where magma undergoes fragmentation. This has led to the utilisation of volatile concentration gradients preserved in melt-filled, open-ended melt inclusions (commonly termed melt embayments or re-entrants) as a geospeedometer, which can record the decompression rates of magmas during their final ascent to the surface (Anderson 1991;Liu et al. 2007;Lloyd et al. 2014;Ferguson et al. 2016;Myers et al. 2018Myers et al. , 2021Georgeais et al. 2021;Saalfeld et al. 2022). Volatile gradients preserved in glassy melt embayments in the crystal cargo of an eruptive deposit reflect the diffusive loss of volatiles within the embayment melt as it strives to achieve and then maintain equilibrium with the surrounding melt phase, which is degassing due to decompression during ascent. It is this disequilibrium that allows the gradients to be used as a geospeedometer for the average decompression rate of the magma in the mid-upper conduit, which can then be related to the magma ascent rate (Liu et al. 2007;Humphreys et al. 2008;Lloyd et al. 2014;Ferguson et al. 2016;Myers et al. 2018Myers et al. , 2021Georgeais et al. 2021;Saalfeld et al. 2022). Information about magma storage conditions from the volatile contents of sealed melt inclusions (e.g. Liu et al. 2006;Johnson et al. 2011;Druitt et al. 2016) can be used as starting conditions in these decompression models, which through reproducing the volatile concentration gradient profiles measured in the melt embayments can provide timescales for magma ascent. Together, these are useful parameters for estimating the warning time a volcano may give prior to a given explosive eruption, with times reflecting the onset of magma movement towards the surface, when sudden changes in the behaviour of the volcano may be expected (e.g. increased seismicity, volcanic tremor, ground deformation, increased gas flux, thermal changes).
So far, studies of magma ascent using this method have largely been confined to explosive mafic (e.g. Lloyd et al. 2014;Ferguson et al. 2016;Newcombe et al. 2020) or large, caldera-forming felsic eruptions (e.g. Myers et al. 2018Myers et al. , 2021Geshi et al. 2021;Saalfeld et al. 2022). Larger rhyolite eruptions (including those up to super-eruption scale, i.e. >450 km 3 DRE) have been shown to have similar magma ascent rates, regardless of the initial pattern of magmatic unrest (e.g. whether the initial explosive activity is continuous or episodic: Myers et al. 2018). Although smaller explosive rhyolitic events (i.e. <20 km 3 DRE) are far more common in the geological record (e.g. Jurado-Chichay and Walker 2001;Castro and Dingwell 2009;Barker et al. 2015), they have so far been neglected in studies of magma ascent dynamics, despite being crucial for future volcano monitoring initiatives. Here, we present the results of an investigation into the pre-eruptive magma storage and the ascent of a series of young (post-25.4 ka) rhyolitic eruptive events from the Ōkataina Volcanic Centre. These results provide a framework of expected decompression rates that can be expected in the region immediately prior to future explosive activity, spanning a range representative of these more common, smaller rhyolitic eruption sizes (from 5 to 17.5 km 3 ). A future eruptive episode at Ōkataina may be preceded by a period of unrest on the order of weeks to years (as the eruptible magma is accumulated and primed to erupt: Elms 2022). However, it is also known that explosive eruptions can occur with very little prior unrest depending on the mechanism of eruption initiation (for example if the trigger does not involve significant magma movement: see Cassidy et al. 2019, for overview). If this was the case, then the magma ascent rates presented here (on the order of 0.5-2 m.s −1 , with final ascent times of tens of minutes to hours), derived from the modelled decompression rates of the magma during ascent through the upper conduit at the onset of the eruption, may provide a lower limit for the amount of time that there might be unrest or geophysical signals prior to a future event

Geological setting
Ōkataina Volcanic Centre ) is a caldera complex situated at the northern end of the rhyolite-dominated Central Taupō Volcanic Zone of Aotearoa New Zealand (e.g. Wilson et al. 2009;Chambefort et al. 2014;Fig. 1a). Magmatism in the Ōkataina area is inferred to have commenced at c. 0.75 Ma, based on zircon U-Pb ages from plutonic fragment-hosted zircon crystals, while the eruptive record begins at c. 0.63 Ma (Manning 1996;Shane et al. 2012). Caldera-forming eruptions occurred at c. 0.56 Ma, 0.32 Ma and 53 ka, alongside numerous intra-caldera events, with the possible further collapse of part of the existing caldera structure at 33 ka Cole et al. 2010Cole et al. , 2014Deering et al. 2011;Shane et al. 2012).
Since the most recent caldera-forming eruption, the c. 53 ka Rotoiti event, magmatic and eruptive activity appears to have been confined to within the caldera boundary and surrounding structural embayments (Jurado-Chichay and Walker 2000;Benson et al. 2021). Post 25.4 ka activity has consisted of nine rhyolitic eruptive episodes and two basaltic eruptions that collectively form the Rotorua Subgroup (Healy et al. 1964;Fig. 1b; Table 1). The rhyolitic episodes generally involved early explosive activity followed by voluminous lava extrusion along two NE-SW-oriented vent lineations, creating Mt. Tarawera and the Haroharo Massif, plus a smaller, NW-SE-oriented dome complex in the Ōkareka Structural Embayment area, southwest of Haroharo (Fig. 1b).
The physical volcanology and geochemistry of these events have been extensively studied (e.g. Nairn 1992Wright 2000;Bowyer 2001;Nairn et al. 2001Nairn et al. , 2004Speed et al. 2002;Kobayashi et al. 2005;Smith et al. 2004Smith et al. , 2005Smith et al. , 2006Darragh et al. 2006;Kilgour and Smith 2008;Rocco 2020;Sas et al. 2021;Elms 2022). However, estimates of the timescales on which the associated pre-and syn-eruptive magmatic processes operated are largely lacking, with current information derived from  (2002), Wilson et al. (2009) andBenson et al. (2021). Eruption ages and total volumes for each eruptive episode are in Table 1 20 Page 4 of 20 numerical conduit modelling or comparison with similar rhyolitic events from other volcanoes (e.g. Papale et al. 1998;Sherburn and Nairn 2004). This study aims to provide additional constraints on pre-eruptive magma storage depths and new insights into the timings and patterns of magma ascent at the onset of these young eruptions by using a geo-speedometer specifically designed to capture these constraints directly from evidence preserved in erupted material (Liu et al. 2007;Myers et al. 2018).

Materials and methodology
The material was sampled from the first-erupted explosive products of eight of the nine Rotorua Subgroup eruptive episodes, with the 25.2 ka Te Rere episode omitted as the basal deposits from this event do not contain quartz. Following the preparation methods of Myers et al. (2018), larger pumice fragments were gently crushed to liberate the crystal cargo, the material was sieved and quartz crystals were picked from the 500-1000-μm fraction. These crystals were then immersed in isopropyl alcohol, so that suitable melt inclusions (Fig. 2a) and melt embayments (Fig. 2b) could be identified and picked. In the 21.9 ka Ōkareka eruption sample, where basalt was co-erupted with the initial rhyolite, care was taken to avoid any inclusions or embayments that visually seem to have been affected (e.g. by partial intrusion or mixing) by the basalt. Suitable inclusions and embayments chosen for analysis were glassy and free of both bubbles and mineral inclusions (e.g. Rose-Koga et al. 2021) and comprised 75-90% of the total inclusions (depending on the eruptive episode) and >90% of the total embayments present in the samples. The presence of bubbles in some inclusions, and a lack of evidence for cracking of the mineral host, indicates that the young Ōkataina magmas were volatile-saturated at the time of quartz crystallisation and melt inclusion entrapment. This indicates that the depth estimates derived from the volatile saturation pressures are reliable if little to no post-entrapment H 2 O modification has occurred. Melt embayments were chosen based on being as straight and cylindrical as possible, ideally with a vapour bubble at the mouth, or, where no such bubble was present, continuity between the melt embayment and selvedge glass (after Liu et al. 2007;Lloyd et al. 2014;Myers et al. 2018;deGraffenried and Shea 2021).
The application of 1D diffusion modelling ideally requires straight, cylindrical embayments as irregular geometries can result in diffusion in multiple directions at once (i.e. 2D or 3D diffusion), and narrowing or widening of the embayment can affect the rate of volatile exchange between the melt embayment and surrounding, decompressing melt (Liu et al. 2007;deGraffenried and Shea 2021). The effects of necking in melt embayments were quantified numerically by deGraffenried and Shea (2021), who showed that there can be significant modification to decompression rates retrieved when necking is greater than 40-50% (i.e. melt embayment neck width/interior width is <0.6), a result we take in consideration when evaluating our embayments.
Melt inclusions (6-14 per eruption) were prepared as double-sided polished wafers and analysed for both H 2 O and CO 2 by transmission Fourier transform infrared (FTIR) spectroscopy. Glass selvedges (representing the host melt or matrix surrounding the quartz crystal) on seven wafers from the Kaharoa (3), Rerewhakaaitu (3) and Ōkareka (1) eruptions were also analysed, to assess the residual volatile contents of the melts at the point of quenching. Measurements were conducted on a Bruker Hyperion 2000 IR microscope connected to a Bruker Vertex 70 spectrometer bench Table 1 Summary of young eruptive episodes at Ōkataina (note that the Te Rere is not included in this study due to the lack of quartz in the earliest explosive deposits, nor are the basaltic Rotokawau and Tarawera events). Magma DRE volumes are for the entire eruptive episode, which in all young Ōkataina rhyolitic events involves both explosive and effusive products, with explosive volumes typically comprising a smaller portion of the total erupted volume. References: 1, ; 2, Walker et al. (1984); 3, Lowe et al. (2013); 4, Beanland and Houghton (1991)  In this equation, c i is the concentration of the relevant volatile species, M i is the mass of the species in g.mol −1 , A is the absorbance (height of the 2350 cm −1 molecular CO 2 peak), ρ is the glass density in g/L (calculated using an iterative approach based on the H 2 O content: Supplementary Table 1), d is the thickness of the wafer in cm, and ɛ is the molar absorption coefficient (L/mol.cm). A ɛ value of 1214 L/mol.cm −1 was used to calculate total CO 2 content from the 2350 cm −1 peak (Behrens et al. 2004). The limit of detection of the CO 2 peak from the background was determined to be c. 10 ppm. Due to oversaturation of the mid-IR total H 2 O peak at ~3570 cm −1 , total H 2 O contents were calculated by combining the near-IR 5230 cm −1 molecular H 2 O and 4520 cm −1 OH − peaks using the formulas and parameters of Zhang et al. (1997). The principles of the Beer-Lambert law  Myers et al. (2018). Negative reflectance FTIR peaks were converted into H 2 O contents after Yasuda (2014) using the equation shown still apply to this method, but the equations of Zhang et al. (1997) account for the molar absorption coefficients of for the 5230 cm −1 molecular H 2 O and 4520 cm −1 OH − peaks being dependent on total H 2 O content, which produces a more accurate total H 2 O content, especially when over-saturation of the main 3570 cm −1 total H 2 O peak occurs.
The total H 2 O and CO 2 contents of melt inclusions were then combined with existing magmatic temperature estimates  to calculate vapour saturation pressures using VolatileCalc (Newman and Lowenstern 2002). Many of the youngest Ōkataina eruptions contain multiple magma types, which have different intensive parameters, such as temperature, which are important for accurate magma ascent modelling. However, not all of these magma types are visually distinguishable. Where this is the case, existing stratigraphic, geochemical and petrological data (e.g. Nairn 1992Nairn et al. 2001Nairn et al. , 2004Speed et al. 2002;Kobayashi et al. 2005;Smith et al. 2004Smith et al. , 2005Smith et al. , 2006Darragh et al. 2006;Shane et al. 2007Shane et al. , 2008Kilgour and Smith 2008) were used to ensure the relevant magma was correctly identified and the appropriate temperature estimates (obtained from the summary of Rotorua Subgroup magmatic properties in Table 1 of Smith et al. 2005, for consistency) were used. Vapour saturation pressures were then converted into approximate magma storage depths using a crustal density of 2640 kg.m −3 , appropriate for the regional Mesozoic metasedimentary basement rocks (Tenzer et al. 2011;McNamara et al. 2014).
Quartz-hosted melt embayments were prone to disintegration when prepared as double-sided wafers, so samples were instead prepared as single-sided, flat polished surfaces and analysed by reflectance FTIR. The same apparatus and set-up was used in reflected light, and 800-1200 scans were measured for both background and samples to reduce noise to acceptable levels. Backgrounds were measured using a gold plate. Although this method precludes the determination of CO 2 contents, the low levels of CO 2 in co-erupted melt inclusions and the frequency of rhyolitic melt embayments modelled in previous studies lacking measurable CO 2 suggest this approach does not introduce significant error (e.g. Myers et al. 2018Myers et al. , 2021Saalfeld et al. 2022). The trough at reflectance ~3650 cm −1 on the reflectance FTIR spectrum (Fig. 3b) was converted to H 2 O concentrations following the method of Yasuda (2014). This study empirically determined a relationship between the depth of the trough normalised to the height (i.e. reflectance value) of the peak baseline and the total H 2 O content, which for rhyolite glass is described by the equation: The concentrations obtained using this method were validated by comparison with the concentrations obtained using the alternative calibrations of Hervig et al. (2003) and Johnson et al. (2011), which relate total H 2 O content to the depth of the trough normalised to the reflectance value of the spectrum at 4000 cm −1 . Both alternative calibrations produce total H 2 O concentrations in close agreement with those produced using the Yasuda (2014) calibration, indicating a lack of instrument dependence to this calibration. Therefore, the Yasuda (2014) calibration was used in this study, as it is the most recent and least complex of the three calibrations. Transects were measured along 2-4 melt embayments per eruption to obtain total H 2 O gradients. These were then combined with the melt inclusion volatile data and previously determined Fe-Ti oxide-derived temperature estimates of 731-762 °C, plus 871 °C for the Rotorua eruption ; Table 2) to create the input parameters required to obtain decompression rates, and associated rates and timescales for magma ascent, in the different eruptions studied. Decompression modelling was carried out following the MATLAB script of Myers et al. (2018), utilising the solubility relationship of Liu et al. (2005) and the diffusivities of Zhang and Behrens (2000), with a model dP/dt (decompression rate) step size of 0.0001.

Results
Analysed sealed melt inclusions were all bubble-free, and round or slightly faceted (e.g. Fig. 2a), ranging from c. 55 to 270 μm in diameter. Sealed melt inclusions (Table 2) from all eruptions (with 2 s.d. uncertainties for each inclusion) have moderate to high total H 2 O contents of 2.92-5.82 ± 0.12-0.40 wt% (Fig. 4a) and low CO 2 contents of 13-126 ± 1-8 ppm (Fig. 4b). The Mamaku and Waiohau samples have the most uniform H 2 O contents, whereas the Rotomā and Rotorua melt inclusions display bimodal H 2 O contents (Fig. 4a). The Rotorua and Mamaku melt inclusions have the most uniform CO 2 values, whereas the most variable are the Whakatāne, Rotomā, Rerewhakaaitu and Ōkareka (Fig. 4b). There is a slight trend towards younger eruptions having lower volatile contents; however, two exceptions are the Waiohau (with relatively low H 2 O) and Rotorua (with relatively low CO 2 ) events. Glass selvedges (Table 3) contain 1.04-1.89 ± 0.11-0.21 wt% H 2 O, with five out of seven selvedges preserving CO 2 contents of 10-16 ± 1-6 ppm, and two returning CO 2 contents below detection (<10 ppm).
VolatileCalc vapour saturation pressures for all measured sealed melt inclusions (Newman and Lowenstern 2002) calculated using these volatile data and the temperature data of Smith et al. (2005) (Table 3). These glass selvedge depth estimates are associated with uncertainties of c. 13-41%, due to the larger relative uncertainties on the volatile measurements. The abundances and shapes of melt embayments in quartz crystals from different eruptions vary widely, and many were unsuitable for this kind of study due to their complex geometries (a common feature in magmatic systems, e.g. Ruefer et al. 2021). Combined with the fragility of the host crystals, the overall number of embayments analysed means insights into individual eruptions are limited, and therefore, our data represent an overview of magma ascent for the explosive rhyolite eruptives of the Rotorua Subgroup as a whole. Melt embayments used for analysis (Table 4) ranged from c. 110 to 340 μm in length. The majority had a vapour bubble at the embayment mouth (e.g. Fig. 2b); however, five (one each from the Whakatāne, Rotomā and Rerewhakaaitu, and two from the Rotorua) lacked any noticeable bubble and instead passed directly into a glassy selvedge. Transects of 5-13 points, of a minimum of a 20×20 μm square spot size (dependent on embayment length, width and the presence of any cracks resulting from fracturing during sample preparation: Fig. 2b), were measured along each embayment from the innermost point to the embayment mouth, with all measured embayments preserving gradients in H 2 O content along their length (Supplementary Table 2). Spectra that indicated any degree of overlap with the host quartz were not used for modelling. Replicate analyses indicate an uncertainty on the H 2 O measurement of c. 0.25 wt%. H 2 O contents preserved in the innermost parts of embayments (i.e. the beginning point of each modelled H 2 O gradient) were often lower than those of the sealed melt inclusions from the same eruption; however, in some cases (e.g. the Whakatāne or Waiohau), they overlapped (Fig. 4a).
Calculation of ascent rates and times from the modelled decompression rates was carried out using a fragmentation pressure of 10 MPa, corresponding to a fragmentation depth of c. 0.5 km, based on the depth estimates from the glass selvedges, assuming that the melt quenched soon after fragmentation at a comparable depth (Table 3).
Owing to the scatter in the H 2 O data for each event, a range of starting H 2 O concentrations and pressures (using the interior embayment H 2 O as an initial guide) was used in modelling, covering the H 2 O ranges of both the sealed inclusions and the inner parts of embayments. Best-fit starting conditions are often derived from initial H 2 O concentrations (and derived pressure estimates) that are slightly higher than those of the inner embayments (discussed below). Best-fit profiles to measured gradients had χ 2 values of 0.10-0.57.
Modelled decompression rates (0.0032-0.044 MPa.s −1 ) for the intra-caldera events (i.e. those erupting from Tarawera and Haroharo) yield calculated ascent rates of 0.10-1.67 m.s −1 , corresponding to ascent timescales of 20 min to 3.9 h (Table 4; Fig. 5). Ascent times and rates from within single eruptions were generally consistent (i.e. within error of each other); however, the two ascent times from the Ōkareka event are significantly different from each other, likely reflecting the differing H 2 O contents preserved in the respective melt embayments. The Rotorua event from the Ōkareka Structural Embayment region (Fig. 5) yielded seemingly more rapid ascent rates of 1.59-4.44 m.s −1 (derived from decompression rates of 0.043-0.117 MPa.s −1 ), with magma rising over the course of 22-34 min. Uncertainties on both ascent rates and ascent times were estimated from the range of dP/ dt (decompression rate in MPa.s −1 ) values that produced profiles with acceptable fits to the measured gradients (i.e. χ 2 <1) and are precise to within a factor of ~1.65 (Table 4; Fig. 5).
Each best-fit model profile was modelled again using the upper and lower bounds of the temperature estimates for each eruption  to estimate the overall likely range of magma ascent rates and timescales represented in each event (presented in the Supplementary  Table 4; Fig. 5). This approach produced a slightly wider (by a factor of ~1.14-1.25) overall range of magma decompression rates (0.003-0.05 MPa.s −1 ), ascent rates (0.08-1.9 m.s −1 ) and ascent times (18 min to 4.7 h) for eruptions from the main intra-caldera region. Likewise, the Rotorua eruption from the Ōkareka Structural Embayment also returned slightly wider ranges of values (by a factor of ~1.17-1.33) when estimating the effect of the temperature uncertainty, with magma decompression rates of 0.04-0.13 MPa.s −1 , ascent rates of 1.5-4.8 m.s −1 and ascent times of 20-36 min. While the magma involved in the Rotorua event appears to have decompressed and ascended more rapidly, two out of three Rotorua decompression and ascent rates overlap with the data obtained for intra-caldera eruptions. As such, whether magma decompression and ascent in the Ōkareka Structural Embayment region is truly more rapid than inside the caldera is equivocal. Perfectly cylindrical melt embayments are rare in the young Ōkataina quartz crystals, and therefore, the majority of the embayments studied here (17 out of 21) exhibit a degree of gradual narrowing towards the embayment   Liu et al. (2007). Black crosses and error bars show the accepted rates and times and the associated uncertainties as presented in Table 4. Coloured bars represent how the range in result values shift (compared to the ranges shown by the crosses) when the upper and lower bounds of the temperature uncertainty ) are taken into account mouth (cf. the sudden "necks" of deGraffenried and Shea 2021). However, four melt embayments had a more pronounced necked geometry with a neck width/interior width ratio of <0.35 (Supplementary Tables 3, 5). The effect of this on the model decompression rates results is discussed below.

Discussion
Determining the rate of magma ascent prior to volcanic eruptions is crucial for volcanic monitoring and for interpreting future unrest episodes. Below, we discuss the trends observed in the melt inclusion data from Ōkataina, the inferred storage conditions prior to eruption and the calculated timescales of ascent as inferred from melt embayments. We then compare our data to other volcanic systems and discuss the implications for eruption mechanisms and volcanic monitoring initiatives at Ōkataina and similar volcanoes worldwide.

Scatter in melt inclusion volatile data
Along with the generalised trend towards lower volatile concentrations and shallower storage depths in younger events seen across the entire dataset, there is also a notable scatter found within individual events (Fig. 6). For H 2 O, this scatter can be interpreted as being caused by a variable combination of two degassing-related processes. The first is that differing volatile contents within sealed melt inclusions could represent entrapment at different times within the process of degassing during magma storage, recharge, crystallisation and slow initial ascent, or within different regions of a system that was heterogenous with respect to volatile contents. The second is the modification of volatile contents via post-entrapment diffusive H loss through the host crystal, where the partial re-equilibration of H 2 O between the inclusion and surrounding, degassing melt is driven by decompression (Fig. 6). This latter process occurs on timescales of days to weeks at typical rhyolitic magma temperatures (e.g. Myers et al. 2016Myers et al. , 2019 and results in misleadingly low H 2 O measurements and subsequent depth estimates, producing a great deal of scatter. This combination of processes could explain the lack of any definitive trends (such as those solely indicative of isobaric degassing) in the volatile and derived pressure data (Fig. 6). Therefore, where final magma storage depth estimates are diverse, the greater pressure estimates are likely to provide the most accurate indication of the true magma starting storage depth.

Pre-eruptive storage conditions at Ōkataina
The range of volatile concentrations measured in melt inclusions in this study (2.92-5.82 wt% H 2 O and 13-126 ppm CO 2 ) agrees well with previous studies of Ōkataina volatiles (e.g. Shane et al. 2007Shane et al. , 2008Johnson et al. 2011) and those of other centres in the TVZ (e.g. Liu et al. 2006;Bégué et al. 2014;Myers et al. 2019) in having moderate to high H 2 O and very low CO 2 (Fig. 7a). The exceptions to this are the CO 2 data of Shane et al. (2007Shane et al. ( , 2008, where melt inclusions in the Rerewhakaaitu and Ōkareka eruptions were found to have up to 1859 and 3359 ppm CO 2 , respectively. These discrepancies may be due to Shane et al. (2007Shane et al. ( , 2008 measuring H 2 O and CO 2 contents in their samples by secondary ion mass spectrometry (SIMS), whereas the above studies  (Newman and Lowenstern 2002) showing scatter indicative of a combination crystallisation, iso-and polybaric degassing and post-entrapment diffusive H loss all used FTIR. If the samples measured in Shane et al. (2007Shane et al. ( , 2008 were mounted in epoxy resin, which can outgas under vacuum during SIMS analysis, this may have resulted in artificially high CO 2 results (Hauri et al. 2002;Esposito et al. 2014;Rose-Koga et al. 2021). The volatile measurements presented here, coupled with relatively cool magmatic temperatures Johnson et al. 2011;Table 2), result in shallow apparent model pressures for pre-eruptive magma storage of  MPa (Newman and Lowenstern 2002), and associated depth estimates of 2.2-7.4 km below the surface (Fig. 4c). Notably, these pressure and depth values overlap with regions in the crust where recent geophysical studies indicate the presence of partial melt beneath the modern volcano . Due to the likelihood of post-entrapment diffusive H loss from the sealed melt inclusions (indicated by the variation in H 2 O: Myers et al. 2016Myers et al. , 2019 as well as either isobaric or polybaric degassing (indicated by the variation in CO 2 alongside the H 2 O), the shallower pressures/ depths measured in each event (which are closely linked to H 2 O content: Table 2) do not accurately reflect the storage depth of the magma, and therefore the highest pressures and depths from each event (derived from the highest-H 2 O sealed melt inclusions) are deemed the most reliable. The resulting accepted magma storage pressures and depths thus range from 117 to 192 MPa and 4.5 to 7.4 km, respectively. Importantly, the overall trend of shallower magma storage depths towards the present day is still apparent in these deeper pressure/depth estimates, trending from 7.3 and 7.4 km in the Ōkareka and Rerewhakaaitu eruptions (respectively) to 5.1 km and 5.0 km in the Whakatāne and Kaharoa events. The Waiohau event, however, does not follow this overall trend, with its sealed inclusions preserving relatively low H 2 O concentrations resulting in an accepted magma storage depth of 4.5 km.

Melt embayments and timescales of ascent
While the limited number of samples (two to four melt embayments for each eruption) precludes detailed analysis of the initial phases of individual eruptive episodes, the data presented here provide new insights into the ascent dynamics of the most recent activity at Ōkataina as a whole, and insight for rhyolitic eruptions overall.
Modelling of H 2 O diffusion to produce best-fit profiles to the measured gradients was attempted using a range of starting H 2 O contents, from both sealed melt inclusions covering the range of values for each eruption (which did not produce acceptable fits), and from the innermost parts of the melt embayments. However, the best fits to the data required initial model H 2 O contents that were slightly higher than those measured in the inner parts of the embayments, by 0.1-1.6 wt% (Table 4). One Ōkareka-episode melt embayment, with an inner measured water content of 2.9 wt%, required a model initial H 2 O content that was 2.6 wt% higher to obtain a best-fit profile to the measured gradient. However, this is the same best-fit model initial H 2 O concentration as was found for the other Ōkareka-episode melt embayment, which has a higher inner melt embayment H 2 O content of 5.0 wt% ( Table 4). The discrepancy between these two Ōkarekaepisode melt embayment H 2 O contents is likely due to the lower-H 2 O embayment being shorter in length (Table 4). As mentioned above, the lack of melt embayment CO 2 measurements in this study is not likely to have a significant effect Fig. 7 Comparison of the H 2 O and CO 2 contents of sealed melt inclusions from the youngest Ōkataina eruptives with data from other rhyolitic eruptives a in the TVZ (Liu et al. 2006;Johnson et al. 2011;Bégué et al. 2014;Myers et al. 2019) and b globally (Newman and Chesner 1989;Wallace et al. 1999;Wallace 2005;Chesner and Luhr 2010;Druitt et al. 2016;Myers et al. 2016Myers et al. , 2019Meszaros et al. 2019;Waelkens et al. 2017Waelkens et al. , 2022. Ōkataina data from Shane et al. (2007Shane et al. ( , 2008, with CO 2 ranging from 13-1859 and 50-3359 ppm, respectively, are omitted due to scale, and data for Chaitén (Castro and Dingwell 2009) are omitted due to having no detectable CO 2 on the model results, as similar studies on other events have shown that rhyolitic melts regularly lack detectable CO 2 , especially in systems where sealed inclusions begin with low values (Myers et al. 2018(Myers et al. , 2021Saalfeld et al. 2022).
As discussed above, sealed melt inclusions in some eruptions show variations in preserved volatile abundances indicative of diffusive H loss as a result of re-equilibration during a slow, initial ascent phase (e.g. Myers et al. 2019). This inference is supported by the melt embayment data, where inner melt embayment H 2 O concentrations frequently overlap with the lower H 2 O measurements of the sealed melt inclusions (e.g. in the Waiohau, Rotorua and Ōkareka data), but are often at the lower end of, or below, the H 2 O ranges seen in the sealed inclusions (Fig. 4a). Additionally, there is a moderate correlation (R 2 = 0.502) between melt embayment length and the difference in H 2 O content between the lowest-H 2 O sealed melt inclusion for an eruption and the inner parts of the associated melt embayments, with shorter melt embayments having the biggest differences ( Fig. 4a; Table 4). These observations imply a two-stage ascent model, with a slower initial phase allowing time at the beginning of the ascent to reequilibrate both melt inclusions and melt embayments to lower H 2 O concentrations, followed by a more rapid, final ascent, the average rate of which is preserved in the melt embayment H 2 O re-equilibration gradients (e.g. Myers et al. 2018Myers et al. , 2019. Shorter embayments appear to have reequilibrated more completely before the more rapid, final ascent began, whereas the inner portions of longer melt embayments preserve H 2 O contents closer to those present in the final stages of magma storage or immediately prior to final magma ascent. These interpretations could be tested further by collecting incompatible trace element data from sealed melt inclusions and the innermost parts of melt embayments to restore the crystallisation-and degassing-modified H 2 O contents to their original abundances at the time of entrapment (e.g. Johnson et al. 2011;Myers et al. 2019).
Ōkataina melt embayments rarely have perfect geometries for 1D diffusion modelling, with 17 out of 21 samples having necking ratios ranging from 0.62 to 1.14, and four samples having more pronounced necking of 0.20-0.32. The degree to which this necking affects modelling results was previously assessed numerically by deGraffenried and Shea (2021), and for natural samples from the Otowi and Tshirege members of the Bandelier tuffs by Saalfeld et al. (2022). Both studies find that the 1D approach is reasonable (within a factor of ~1.8) for embayment geometries with necking ratios >0.6. For instance, Saalfeld et al. (2022) found that even for their most extremely necked sample (44 % necked, i.e. a ratio of 0.56), the 2D model decompression rate was closely aligned with the 1D result, where errors in determined decompression rates were dominated by the large temperature uncertainty (~130 °C, factor of 3-4 offset in extracted decompression rates). However, deGraffenried and Shea (2021) found that differences between 1D and 2D results were offset by a factor of 3-4 when the necking ratio is <0.4. For the 17 embayments with necking ratios of 0.62-1.14, a comparison with the results present in Table 1 of deGraffenried and Shea (2021) suggests appropriate "decompression rate multipliers", or correction factors to account for non-ideal geometry, range from 0.95 to 1.65 (Supplementary Table 5). It is worth noting that the experimental results of deGraffenried and Shea (2021) are not an exact match to the embayments and magmatic conditions studied here, due to the different parameters utilised in the experiments (i.e. a temperature of 800 °C and embayment length of 350 μm, versus the lower temperatures and variable embayments lengths studied here). However, especially when compared against the study of Saalfeld et al. (2022), we believe that the results of deGraffenried and Shea (2021) provide a close approximation to test our embayment geometries against and infer that the effect of geometry in these examples is comparable to the uncertainties derived from both the range of model profiles with acceptable misfits (χ2 <1, Supplementary Table 3) and from temperature uncertainties (Supplementary Table 4). As such, our interpretations are not significantly affected. For the four embayments with more pronounced necks, with ratios of 0.20-0.32, the DRMs are likely to be much higher. Estimated correction factors for these four samples range from 8.70 to 14.80, despite the 1D extracted model rates presented here being consistent with other profiles from embayments within the same eruption (Table 4 Table 3). While attempts at 2D decompression modelling to fully quantify the geometry effect were beyond the scope of this study, these estimates combined with estimates of the impact of the temperature uncertainties indicate that the model decompression rates (and ascent rates and times) are reliable to within a factor of 1.5-2.
All the intra-caldera events studied here have broadly similar final magma ascent rates and timescales. There are no notable temporal trends in the model results, nor are there any marked differences between the rates obtained for eruptions from Tarawera (0.006-0.036 MPa.s −1 , 0.19-1.35 m.s −1 , 20-230 min) vs. Haroharo (0.0032-0.044 MPa.s −1 , 0.10-1.67 m.s −1 , 23-140 min). This indicates that similar ascent dynamics occurred on both sides of the intra-caldera magmatic system. However, magma ascent during the Rotorua event from the structural Ōkareka Structural Embayment region may have occurred more rapidly, with decompression and calculated ascent rates ranging from being comparable with the intra-caldera events up to 0.12 MPa.s −1 and 4.4 m.s −1 , respectively ( Fig. 5; Table 4). This may reflect the fact that the initial Rotorua magma was considerably hotter than the other events, at 871 °C vs. 731-762 °C for the intra-caldera events , meaning that the magma would be less viscous and able to ascend more rapidly.
While the timescales of ascent for the Rerewhakaaitu, Rotorua, Waiohau, Rotomā, Mamaku, Whakatāne and Kaharoa events are all within uncertainty, one timescale determination for the Ōkareka event is significantly longer (Fig. 5; Table 4). The reasons for this are unclear, given that the melt embayment with the profile producing the longer timescale is the shorter of the two, and both require the same initial conditions to achieve the best-fit model output. Interestingly, however, ascent rates obtained from both Ōkareka-episode melt embayments are no different than those for the other intra-caldera events, despite the initial Ōkareka rhyolite being co-erupted with (hotter, less viscous) basaltic material that contains high-Mg mantle-derived olivine (Barker et al. 2020). This indicates that the presence of subordinate volumes of basaltic magma in the initial stages of this rhyolitic event did not exert a major control on the rate of rhyolitic magma ascent.

Comparisons with other rhyolitic systems
The volatile contents of the sealed melt inclusions analysed in this study are typical of rhyolitic eruptives from the Taupō Volcanic Zone (Fig. 7a). The range of H 2 O contents measured is also comparable to that of many other rhyolitic centres worldwide; however, the abundance of CO 2 is relatively low (see also Johnson et al. 2011;Fig. 7b).
Many of the ascent rates determined in this study are deemed relatively slow for explosive rhyolitic events by some authors (e.g. Rutherford 2008;Gonnermann and Manga 2013). However, as noted by Myers et al. (2018), values in these review papers are mostly derived from numerical conduit models (e.g. the 5-8 m.s −1 of Papale et al. 1998), which represent the final and rapid ascent rates leading up to fragmentation. In contrast, volatile diffusion modelling of melt embayment gradients (Liu et al. 2007;Lloyd et al. 2014;Ferguson et al. 2016;Myers et al. 2016Myers et al. , 2018Myers et al. , 2021Georgeais et al. 2021;Saalfeld et al. 2022) is specifically designed as a geospeedometer for average magma ascent from source to surface, and so rates derived from such models are likely to be more representative of the total magma rise process.
Comparing the decompression and ascent rates obtained in this study with those modelled for other rhyolitic events, or estimated from direct observations, presents some important implications for rhyolitic eruption dynamics. The only directly observed rhyolitic event of a similar size to the initial explosive events of the young Ōkataina eruptions is the May 2008 eruption of Volcán Chaitén in Chile, where physical indicators of magma movement, such as earthquakes, were only felt 24 h before the first ashfall (Major and Lara 2013). This short warning period was taken to indicate that magma ascent was fairly rapid, and this inference is supported by petrological and experimental data suggesting that the magma ascended from >5 km depth to the surface in the space of 4 h, at rates of 0.5-1 m.s −1 (Castro and Dingwell 2009). The ascent time of about 4 h is comparable to the longer magma ascent times obtained for the young Ōkataina events, with 10 out of 21 melt embayment gradients producing model ascent times of 1.0-3.9 h. The estimate of >5 km storage depth of the Chaitén magma is also comparable to the 4.5-7.4 km magma storage depths obtained here for Ōkataina (Tables 2, 4), which in turn are consistent with the depths inferred for the current magma reservoir below Haroharo .
While the Chaitén eruption was similar in eruptive style to the early activity of the young Ōkataina eruptive episodes, ascent rates obtained for super-sized (>450 km 3 magma) eruptions are, interestingly, also very similar (Fig. 8). Model ascent rates obtained for the Bishop Tuff (Myers et al. 2018;Jollands et al. 2020), Huckleberry Ridge Tuff and Ōruanui Formation (Myers et al. 2018), the Bronze Age (Minoan) eruption from Santorini (Myers et al. 2021) and Otowi and Tshirege members of the Bandelier Tuff (Saalfeld et al. 2022) all fall within broadly the same range. The majority of data from these very large eruptions represent ascent rates of <2.5 m.s −1 , similar to the intra-caldera events from Ōkataina and the Chaitén eruptions, with more scattered examples up to >5 m.s −1 , similar to the Rotorua event (Fig. 8). Given the similar depths of long-term magma storage, ascent times in these larger events appear to be similar to the smaller Ōkataina and Chaitén events (Fig. 8), with final magma ascent during the initial phases of the Bishop, Ōruanui and Huckleberry Ridge eruptions occurring over <1-2.5 h (Myers et al. 2018). The model results for the Bandelier Tuff, where a good fit to the data (i.e. χ 2 <1) was achieved, also produced comparable average decompression and ascent rates of 0.041 MPa.s −1 and 1.63 m.s −1 (Otowi Member) and 0.018 MPa.s −1 and 0.47 m.s −1 (Tshirege Member). These decompression and ascent rates correspond to model timescales of a few tens of minutes up to c. 8.5 h (Saalfeld et al. 2022).

Implications for eruption mechanisms and volcanic monitoring
Explosive rhyolitic events range widely in size and vary greatly in the patterns of their initial behaviour, such as the degree of episodicity in their initial explosions, and whether they stall between their slow initial ascent and final rise to the surface (e.g. Wilson and Hildreth 1997;Wilson 2001;Myers et al. 2016Myers et al. , 2018Swallow et al. 2019). However, comparisons between the youngest Ōkataina eruptions studied here and data published for other events of widely varying eruptive volume show that the volume of magma being erupted, and the overall pattern of magma ascent, has practically no bearing on the rates or timings of magma rise through the mid-upper part of the conduit prior to fragmentation (Fig. 8). These similarities are seen despite the potentially large differences in rates and timescales of the initial, slower ascent phase between different events, and differences in overall eruptive style. Likewise, the presence or absence of a mafic magma trigger does not appear to have a notable effect on the rhyolitic ascent rate, as demonstrated by comparisons of the Ōkareka and other Tarawera-vented eruptive episodes vs. those from Haroharo, which show no direct evidence for basaltic input, or between the Ōruanui eruption vs. the Bishop Tuff (e.g. Hildreth and Wilson 2007;Rooyakkers et al. 2018). It is also common for decompression rates obtained from the same stratigraphic unit to be variable by more than an order of magnitude (Figs. 5, 8;Saalfeld et al. 2022). These observations overall imply that the rate of magma ascent recorded in melt embayment-hosted volatile gradients is controlled by factors other than magma volume or overpressure. These factors may include variable conduit conditions and the region of the conduit in which the host crystal ascended (i.e. the middle or the side), or (to a lesser extent: deGraffenried and Shea 2021) the degree to which the external melt around each crystal or near the mouth of the embayment specifically follows an equilibrium degassing path as the magma rises. The difference in ascent rate between the hotter, low-crystallinity Rotorua magma and the cooler, often higher-crystallinity intra-caldera eruptives also suggests that magma viscosity may play a role in ascent rate. However, many eruptions (such as the large events cited here) record larger magma ascent ranges from within single compositional and stratigraphic units, which suggests that conduit dynamics and external controls may override the physical properties of the magma to a certain extent.
The similarity in ascent rates between eruptions encompassing almost the entire recorded size range of rhyolitic explosive events also has important implications for monitoring efforts. Both the magma ascent rates modelled and eyewitness accounts of pre-eruptive tremors during the Chaitén event (Castro and Dingwell 2009) show that rhyolitic eruptions can have rapid onsets with a little warning once magma starts migrating towards the surface. While large or super-sized rhyolite events would be likely to be preceded by warning signs such as higher levels of seismic activity and ground deformation, for small to moderate events, little prior warning may be given (e.g. Ellis et al. 2022). Ōkataina Volcanic Centre (monitored by GNS Science/GeoNet) regularly experiences swarms of earthquakes that have been variably interpreted as fault movement, aqueous fluid migration and/or dike emplacement (e.g. Benson et al. 2021;Hamling et al. 2022). Like other caldera systems, one of the challenges in monitoring Ōkataina is deciphering the context of activity given the setting of the volcano in an active rift (e.g. Berryman et al. 2022;Villamor et al. 2022), versus activity related to the volcano that may signify volcanic unrest. From our analyses of past eruptions, it is clear that pre-eruptive processes likely involve the accumulation of melt at shallow levels (c. 4-8 km depth) prior to relatively rapid ascent to the surface via diking (on the order of a few hours). Magma accumulation at this depth may be difficult to detect, depending on the location of emplacement, depth and volume of magma (e.g. Ellis et al. 2022). In addition, seismicity within the magma system may be limited due to the ductile nature of the magma reservoir (e.g. Illsley-Kemp et al. 2021). Diking of magma to shallower levels would be expected to produce more obvious signals that could be detected by current monitoring, including shallowing seismicity, volcanic tremor, ground deformation and increasing gas flux (e.g. Sparks 2003;Acocella et al. 2015). However, given the rapid ascent rates of past eruptions we infer from our analysis, it may be challenging to interpret and effectively communicate monitoring data prior to the onset of a future explosive eruption. For example, short-term ground deformation signals may be missed by routine InSAR analysis, with a temporal resolution of days to weeks across this region. Combined with a relatively sparse GNSS network, this means that pre-eruptive ground deformation in the Ōkataina region may be difficult to detect (Benson et al. 2021). Therefore, only limited warning times may precede future small-to moderate-volume eruptions from Ōkataina, highlighting the need for detailed analysis of unrest signals and rapid response by monitoring agencies.