Connecting the Dots: the Lava Domes' Perspective of Magmatism Related to an Ignimbrite Flare-Up

Abstract

Dome eruptions associated with rhyolitic calderas offer an important insight into how extremely large (> > 10 km³), rhyolitic magma systems are constructed through time. We focus on rhyolitic calderas in the central Taupō Volcanic Zone leading to, during, and immediately following the 350- to 240-ka ignimbrite flare-up. We identified 103 dome eruptions that are dated between ca. 650 and 150 ka and collated 239 literature whole-rock compositions from these domes. For each composition, we modelled the pressure of magma extraction from the magma mush and the mineral assemblage of the mush using the rhyolite-MELTS geobarometer. We calculated extraction temperatures using zircon saturation geothermometry. We show that magmas are extracted from typically quartz-bearing magma mush at a wide range of depths (~50–425 MPa, ~2–16 km) and temperatures (~750°C to ~850°C). Throughout the central TVZ, there are two dominant extraction pressure modes at 1) 150–175 MPa and 2) 250–325 MPa, consistent with 1) the depth of the brittle–ductile transition (~6 km) and just below typical pre-eruptive storage depths of other TVZ magmas (100–150 MPa, ~4–6 km); and 2) partial melt regions imaged below ~8 km by previous geophysical studies. In some regions, there is a clear correlation between crustal structures, the depth of magma extraction, and the composition of the magmas. In the Whakamaru caldera, the domes erupted inside the caldera following caldera collapse are extracted from ~225 to ~350 MPa at ~810°C and have orthopyroxene-bearing compositions dissimilar to the caldera-forming eruption. These domes are aligned along normal faults, suggesting that rifting creates pathways for magma extraction from a deeper mush rejuvenated by recharge. The domes erupted along the structural margins of the Whakamaru caldera have very evolved, hornblende-bearing compositions, similar to the caldera-forming eruption and shallow, colder extraction from ~100 to ~200 MPa at ~770°C, suggesting the mush feeding these domes is a remnant of the older caldera-forming magma system mobilized along the caldera-bounding faults. Two structural levels of magma extraction at ~6 km and 9 to 12 km are persistent throughout the flare-up period and across the central TVZ region, demonstrating the need for further investigation into the factors controlling the depth of mush development.

INTRODUCTION

Numerous studies have highlighted the complicated, vertically extensive magma plumbing systems that feed the largest (> VEI 7), silicic, caldera-forming eruptions (Dunbar & Hervig, 1992; Annen et al., 2005; Cooper et al., 2012; Jaxybulatov et al., 2014; Cashman et al., 2017; van Zalinge et al., 2017; Cruden & Weinberg, 2018; Kennedy et al., 2018; Black & Andrews, 2020; Boro et al., 2020; Sparks et al., 2021). However, the smaller eruptions from these same caldera systems give us essential insight into the growth of these large magma systems (the ‘caldera cycle’ e.g. Bouvet de Maisonneuve et al., 2021). The existence of these small eruptions raises the questions: why do magmas sometimes assemble to form very large (>100 km³) eruptible bodies, and why do they sometimes erupt as small (<10 km³) eruptions? How are the small and large eruptions related to each other and what do these small eruptions tell us about the state of the caldera system?

In the build-up to caldera-forming eruption, the composition of dome eruptions may record the development of a large, persistent, compositionally mature and homogenous magma system (Hildreth, 2004, 2021; Simon et al., 2007; Troch et al., 2017). Pre-caldera dome eruptions may instead record independent pulses of volcanism largely unrelated to the main caldera-forming events. After the caldera-forming eruption, dome eruptions may represent residual magma left behind after the caldera-forming event (Mucek et al., 2017; Chesner et al., 2020; Tavazzani et al., 2020; Mucek et al., 2021). Alternatively, post-caldera dome eruptions may record intrusion of new mafic magma into the caldera magmatic system (Bachmann et al., 2012; Gelman et al., 2013; Wilcock et al., 2013; Barker et al., 2015; Hildreth et al., 2017). The structural reorganization of the crust during caldera formation may promote pathways for mafic recharge to establish new fractionation and assimilation zones, in which case post-caldera domes may be genetically unrelated to the caldera-forming magma systems (Spell et al., 1993; Girard & Stix, 2009; Pritchard & Larson, 2012; Wilcock et al., 2013). Post-caldera formation, magmatic systems may be much more spatially heterogeneous than the large-scale caldera-forming magma systems that preceded them, with multiple sub-systems developing (Spell et al., 1993; Wilcock et al., 2013). Understanding the lifecycle of caldera-forming systems enables better interpretation of monitoring signals, better assessment of caldera unrest, and ultimately better eruption forecasting (Lowenstern et al., 2006; Kennedy et al., 2018).

In this work, we model the pressure and temperature of magma extraction from the mush Gualda et al. (2019b) to investigate the magma plumbing system of small dome-forming eruptions in the Taupō Volcanic Zone (TVZ), New Zealand. We focus on effusive volcanism leading to, during, and immediately following one of the TVZ's most active periods—the 100-kyr category 3 ignimbrite flare-up (see Gravley et al., 2016). We compare the extraction pressure we obtain for these eruptions to those from large caldera-forming eruptions and consider the magmatic cycle of the caldera plumbing systems, focussing on the depth of the magma mush that fed the dome eruptions across space and time. We emphasize that pre-eruptive storage pressures (Gualda & Ghiorso, 2014; Gualda et al., 2018; Smithies et al., 2023) cannot be obtained for the dome eruptions, due to the extensive syn-eruptive crystallization; we, instead, focus on the depths of extraction, which reveal the distribution of magma mush from which dome magmas were extracted. The dome eruptions give us much better spatial resolution than the caldera-forming eruptions, enabling us to explore spatial relationships between magma extraction depth and volcanic–tectonic features. We compare our results with existing geophysical and monitoring data for the current day TVZ, as well as with models based on petrologic data for the past TVZ.

Definitions

In this study, we define ‘magma mush’ as relatively crystal-rich magma with a rigid framework of crystals and interstitial melt, possibly with bubbles of exsolved volatiles; while magma mush is not typically mobile, it can behave as a fluid under high strain rates (Figure 1) (Bachmann & Bergantz, 2004; Hildreth, 2004; Bachmann & Bergantz, 2008; Cashman et al., 2017). ‘Mush’ in our usage is a description of the physical state of the magma and does not denote a petrogenetic origin (i.e. partial melting vs. fractional crystallization). ‘Eruptible magma’ is defined here as magma that has no finite rigidity and that is typically able to flow like a fluid on eruptive timescales; it is typically melt-dominated magma, possibly with crystals and bubbles. Eruptible magma is extracted (physically separated) from the magma mush prior to pre-eruptive storage of the eruptible magma Gualda et al. (2019b); ‘extraction’ is equivalent to the ‘launching point’ of Blundy (2022). We define pre-eruptive storage as the final accumulation of the eruptible magma prior to eruption. Pre-eruptive storage is represented petrologically by the final equilibration of melt and phenocrysts (Figure 1). ‘Magma bodies’ are discrete parcels of magma (either magma mush or eruptible magma) that are surrounded by country rock or solidified plutonic rocks. The ‘magmatic system’ comprises the mush and eruptible magma involved in a volcanic event, and may include multiple magma bodies.

TAUPŌ VOLCANIC ZONE GEOLOGICAL BACKGROUND

The Taupō Volcanic Zone (TVZ) is a northeast-trending volcanic arc in the North Island of New Zealand (Figure 2). In the last 1.6 Ma, the central portion of the TVZ has produced at least 25 caldera-forming eruptions ranging in eruptive volume from 10s of km³ to super-eruptions of >1000 km³ (Wilson, 1985; Wilson et al., 1986; Nairn et al., 1994; Kampt et al., 1995; Wilson et al., 2006; Gravley et al., 2007; Wilson et al., 2009). In the central TVZ, small to medium volume (0.01–15 km³) explosive and effusive eruptions are more than thirty times as frequent as the caldera-forming eruptions (Kósik et al., 2020). These small eruptions are dominated by rhyolite, particularly effusive rhyolitic lava domes (Figure 2) (Brown, 1994; Wilson et al., 1995; Wilson et al., 2009; Leonard et al., 2010).

Voluminous rhyolitic magmatism in the central TVZ is driven by a complex interplay of arc volcanism related to the subduction of the Pacific Plate under the Australian Plate (Wilson et al., 2009; Reyners, 2013), as well as high rates of tectonic extension in the Australian Plate (~10 mm/yr; Spinks et al., 2005; Wallace et al., 2004). Volcanism and tectonic rifting have a strong spatial and temporal correlation (Berryman et al., 2008; Villamor et al., 2011; Berryman et al., 2022; Villamor et al., 2022), the vents of both small effusive eruptions and caldera-forming eruptions are largely defined and controlled by the tectonic structures (Beresford & Cole, 2000; Spinks et al., 2005; Gravley et al., 2007; Downs et al., 2014b; Kósik et al., 2020).

Crustal structure

Magma systems of the central TVZ span much of the crust. Due to tectonic extension, the crust in the central TVZ is anomalously thin compared to the surrounding crust (Reyners et al., 2006; Wilson & Rowland, 2016). The seismic velocity structure of the central TVZ suggests that the crust is quartzo-feldspathic to a depth of 15 to 20 km (Harrison & White, 2004; Stratford & Stern, 2004, 2006; Harrison & White, 2006). A velocity increase at 15 to 20 km is interpreted either as the elevated Moho below which is mafic underplating and the mantle (Stern et al., 2006), or as a boundary with mafic, heavily intruded and underplated lower crust (Harrison & White, 2004; Harrison & White, 2006). The transition between the seismogenic upper crust and aseismogenic mid to lower crust in the central TVZ is at approximately 6 km depth (Bryan et al., 1999). Below this transition, magnetotelluric models reveal zones of partial melt (<4%) below 8–10 km across the central TVZ. Highly conductive zones are currently observed under the Whakamaru and Ōkataina calderas, suggesting the presence of up to ~50% partial melt (Ogawa et al., 1999; Heise et al., 2007; Heise et al., 2010; Walter, 2014; Heise et al., 2016). Low shear-wave velocities at ~6 to 16 km depth, indicating zones of partial melt near the Rotorua and Reporoa calderas, are in general agreement with these magnetotelluric studies (Bannister et al., 2004). Mush may develop even shallower—plutonic lithics erupted in the Ōkataina caldera complex have compositions and textures that suggest they are crystal cumulates developed at ~150 MPa (~6 km assuming a crustal density of 2.7 g∙cm⁻³; Stagpoole et al., 2020) (Graeter et al., 2015). Pre-eruptive magma storage in the central TVZ may either be contiguous with the mush (Brown et al., 1998b; Deering et al., 2011b; Gualda et al., 2019b; Barker et al., 2020b; Smithies et al., 2023) or the eruptible magma may be extracted from the mush and moved to a separate, shallower pre-eruptive storage region as pods, sills or laccoliths of eruptible magma (Cole et al., 2014; Gualda et al., 2019b; Pamukçu et al., 2021; Smithies et al., 2023). Petrologic data indicates that central TVZ rhyolites are typically stored at ~100–150 MPa (4–6 km), but pre-eruptive storage can be as deep as 11 km and as shallow as 2 km (Bégué et al., 2014b and references therein).

Rhyolite compositional endmembers

Rhyolite compositions in the central TVZ fall along a continuum of compositions. The mineralogy of the rhyolites ranges from relatively crystal-rich (up to 40 wt %), amphibole ± biotite-bearing to crystal-poor (<10 wt %), orthopyroxene-bearing (Ewart, 1967). These mineral assemblages correlate with variations in whole-rock major and trace-element compositions (Deering et al., 2008). Deering et al. (2010) defined the crystal-rich, amphibole ± biotite-bearing endmember ‘R1’; which contrasts with the crystal-poor, orthopyroxene-bearing endmember ‘R2’; some compositions fall in-between the two endmembers, referred to as ‘R1 + 2’. In this study, we use this nomenclature to classify and compare rhyolite compositions (Table 1). We emphasize that ‘crystal-rich’ magmas (20–30 wt %) in the TVZ context are still within the range of crystal contents observed by crystal-poor systems worldwide, and the TVZ magmas do not share the attributes of well-known crystal-rich systems (e.g. Fish Canyon Tuff, total crystal content >40%; Bachmann et al., 2002). In general, the R1 rhyolites are more oxidizing (f_O2 = ΔQFM +1.5 ± 0.3 log units) whilst the R2 rhyolites are more reducing (f_O2 = ΔQFM +0.5 ± 0.3 log units) (calculated by Deering et al. (2010) using 211 Fe–Ti oxide pair compositions and the method of Ghiorso & Sack (1991)). The presence of amphibole in R1, as well as melt inclusion volatile contents, shows that the R1 source material was relatively more hydrous than R2 (Bégué et al., 2015). It is debated whether the compositional variability is driven by fractional crystallization and assimilation processes (Deering et al., 2008; Deering et al., 2010), or whether they are driven by the subduction interface and the mantle source (Deering et al., 2010; Deering et al., 2011a; Rooney & Deering, 2014). With our new geobarometry data, we can explore the crustal controls on these magma compositions by showing the relationship between the magma endmembers and the depth of the final fractionation in the magma mush.

Ignimbrite flare-up

Between ca. 350 and 240 ka, the central TVZ had an ‘ignimbrite flare-up’, a period of above-background numbers of caldera-forming eruptions (Table 2) (Gravley et al., 2016). The ignimbrite flare-up began with the eruption of >2200 km³ magma (dense rock equivalent, DRE) in the Whakamaru Group ignimbrites between 350 and 340 ka (Downs et al., 2014b). Zircon from Whakamaru pumice dated to ca. 600 ka indicates that magmatism in the Whakamaru region took place as early as >250 kyr before the ignimbrite eruptions (Brown & Fletcher, 1999). Gravley et al. (2016) treat this pre-Whakamaru magmatism as the beginning of a ‘magmatic flare-up’ that culminated in the assembly and eruption of the Whakamaru magma system. Following the Whakamaru caldera-forming eruptions, there was voluminous extrusion of rhyolitic lava including the intra-caldera Maroa dome complex and the western and north-western dome complexes (WDC and NWDC, respectively) on the caldera margins (Figure 2b). The volume of these lavas is significant, Leonard (2003) calculated a total volume of ~600 km³ (DRE) for the caldera infill lavas and associated pyroclastic deposits erupted in the first ~60 kyr post-caldera formation.

In the same ~100 kyr period as the Maroa domes extruded, six further caldera-forming eruptions occurred from five caldera centres summarized in Table 2 (Matahina, Chimp, Kaingaroa, Pokai, Ohakuri, and Mamaku ignimbrites) (Gravley et al., 2016). These were an order of magnitude smaller than the Whakamaru eruption; each eruption has estimated magma volumes of 50 to 150 km³. Both the domes and the post-Whakamaru caldera-forming eruptions record a compositional shift in the regional magmatic system. The Whakamaru Group has typical R1 compositions, but the subsequent eruptions move towards R2 compositions with time.

Compared to the Whakamaru caldera, there are relatively few domes preserved at the surface of the other six calderas (Figure 2b). Their relative absence is partially due to burial by sedimentation enhanced by active rifting, burial by subsequent eruptions, and destruction by younger caldera-forming eruptions (Cole et al., 2010; Cole et al., 2014; Downs et al., 2014a; Rosenberg et al., 2020).

In this study, we focus on the dome-forming eruptions in the build-up to, during, and immediately after the 350- to 240-ka ignimbrite flare-up. Whilst the magmatic systems of the caldera-forming eruptions have been well studied (Brown et al., 1998b; Beresford & Cole, 2000; Deering et al., 2011b; Bégué et al., 2014a; Gualda et al., 2018; Gualda et al., 2019b; Smithies et al., 2023), relatively few have focused on the magmatic systems feeding the dome-forming eruptions (Bowyer, 2002; Leonard, 2003; Ashwell et al., 2013). The abundance of dome eruptions associated with the ignimbrite flare-up period—especially when compared to only eight caldera-forming eruptions—gives us a detailed picture of the magmatic systems through time and space. We examine over one hundred dome eruptions and assess structure and evolution of magma systems across the entire central TVZ, with the goal of assessing the temporal evolution of the magmatic systems between the major caldera-forming eruptions. We use dome vent locations to explore the relationships between the subsurface plumbing systems, caldera structures, and tectonic structures. With this new temporal and spatial resolution, we refine our understanding of the magmatic and tectonic processes driving ignimbrite flare-ups.

METHODS

Dome compositions

We compiled a database of whole-rock major-element compositions for dome samples in the central TVZ from published sources and unpublished theses (Supplementary Material A). All samples (n = 239) were analysed by x-ray fluorescence spectrometry for major- and trace-element chemistry (XRF). Published data sources are Deering et al. (2008); Downs et al. (2014a). Unpublished sources with no available methodology are Cole (1979) via GNS Science (2004) analysed at Victoria University of Wellington (four compositions); and Challis (1971) via GNS Science (2004) analysed at an unknown laboratory (one composition).

The remaining data are collated from unpublished theses (Brown, 1994; Beresford, 1997; Richnow, 1999; Milner, 2001; Bowyer, 2002 ; Leonard, 2003). All whole-rock compositions were analysed at the University of Canterbury. Sample preparation was by the individual authors. In all studies rock samples were prepared by washing, drying, and milling in a tungsten-carbide ring mill. XRF analyses were all performed by S. Brown, following the method given by Norrish & Hutton (1969) with modification after Harvey et al. (1973); Norrish & Chappell (1977); Schroeder et al. (1980). The same set of 35 (unspecified) international standards were used for calibration. XRF analyses for Brown (1994); Richnow (1999) were collected on a Philips PW 1400 automatic X-ray fluorescence spectrometer; Beresford (1997); Bowyer (2002); Leonard (2003); Milner (2001) were collected on a Philips PW 2400 automatic X-ray fluorescence spectrometer. Additional trace-element and rare-earth element data were collected by Bowyer (2002) on a subset of samples; these were determined by laser ablation-inductively coupled plasma mass spectrometry (LA-ICPMS) at the Research School of Earth Sciences, Australian National University, Canberra.

We selected rhyolite (>70 wt % SiO₂) domes erupted from ca. 650 ka (before the ignimbrite flare-up, in the period Whakamaru magmatism was active), to ca. 150 ka (100 kyrs after the last flare-up caldera-forming eruption) (Figure 2b, Supplementary Material B & C). We have not included dacites (<70 wt % SiO₂) as they are rare in the central TVZ (Reid & Cole, 1983); only two dacite domes erupted between 650 and 150 ka (Maungakakaramea and Maungaongaonga). Additionally, some TVZ dacites are shown to originate from magma mixing processes (Cole et al., 2001; Millet et al., 2014) and thus do not meet our assumption of the whole-rock composition representing the extracted melt composition. Many of the domes in the central TVZ have not been dated; in these cases, we use the mapped ages based on stratigraphic relationships in Leonard et al. (2010) (Supplementary Material B). A total of 239 compositions from 103 domes were compiled, with sample locations shown in Figure 2b.

We classify each dome into a rhyolite magma endmember (R1, R1 + 2, and R2) based on the criteria of Deering et al. (2008); Deering et al. (2010): firstly by crystal content and mineralogy as reported in the literature, then by FeO_T wt %/MgO wt %, Sr, Y, and Zr concentrations (Tables 1 & 2). Each dome is assigned to the caldera it is associated with, based on location and literature (Bowyer, 2002; Leonard, 2003). Of the 239 compositions and 103 domes, 148 of the compositions come from 54 domes within and on the boundary of the inferred Whakamaru caldera. Because of the abundance of these data, we further subdivide Whakamaru domes into the Maroa volcanic complex, the north-western dome complex (NWDC), and the western dome complex (WDC) (Figure 2b, Supplementary Material B) (Leonard, 2003).

Rhyolite-MELTS geobarometry

Understanding the depth and extent of magma mush is an ongoing challenge in volcanic petrology (Ellis et al., 2014; Graeter et al., 2015; Rubin et al., 2016; Bertolett et al., 2019; Foley et al., 2020; Blundy, 2022). Gualda et al. (2019b) introduced a new method for finding the pressure (correlated to depth) of magma extraction from the mush using the rhyolite-MELTS geobarometer. In their approach, illustrated in Figure 1, Gualda et al. (2019b) assume that the whole-rock composition of the erupted magma represents the composition of the melt that was extracted from the mush. Using rhyolite-MELTS, they find the pressure at which this whole-rock melt composition was in equilibrium with hypothetical mush mineral assemblages. We refer to this as the ‘extraction pressure’, which can be used to infer mush depth. The extraction geobarometer differs from previous applications of the rhyolite-MELTS geobarometer (Gualda & Ghiorso, 2014; Bégué et al., 2014b; Pamukçu et al., 2015; Gualda et al., 2018; Harmon et al., 2018; Pamukçu et al., 2020), which find the ‘pre-eruptive storage pressure’ of magmas by searching for the pressure at which the glass composition is in equilibrium with the phenocryst mineral assemblage coexisting with glass in the studied samples. Pre-eruptive storage pressures from rhyolite-MELTS geobarometry compare well (<100 MPa difference) to results from amphibole geobarometry and H₂O–CO₂ solubility models on the same eruptions (Bégué et al., 2014b; Pamukçu et al., 2015; Gualda et al., 2019a; Pamukçu et al., 2021). We note that these geobarometers measure the pressures of different processes (i.e. quartz + feldspar equilibration with the melt; amphibole crystallization; and volatile saturation). In effusive eruptions, we cannot assume that the glass composition represents the stored melt due to the high probability of syn-eruptive crystallization; we thus cannot calculate pre-eruptive storage pressures, even if glass is preserved in the studied samples. However, even in the case of effusive eruptions, the whole-rock composition of the lava is still representative of the extracted magma, and rhyolite-MELTS extracted pressures can be calculated following the same procedures as those in Gualda et al. (2019b).

The geobarometer uses the rhyolite-MELTS thermodynamic engine (rhyolite-MELTS_v1.0.x; Gualda et al., 2012) with the MELTS-Excel interface (Gualda & Ghiorso, 2015) to calculate the temperature of mineral saturation at a range of pressures (Gualda & Ghiorso, 2014). We use the 239 whole-rock dome compositions as the extracted melt compositions. We run the calculation for each composition from 1100°C to 700°C at 1°C steps and from 500 to 25 MPa at 25 MPa steps. Oxygen fugacity (f_O2) affects the saturation of ferromagnesium phases (Harmon et al., 2018), so we run each composition under a variety of f_O2 values, fixed relative to the QFM buffer, from QFM to ΔQFM+2 in 0.5 log intervals. This encompasses the range of f_O2 values typical in TVZ rhyolites (e.g. Brown et al., 1998b; Smith et al., 2005; Deering et al., 2010; Ghiorso & Gualda, 2013; Barker et al., 2015). Due to uncertainties on the volatile content and, particularly, on the water activity relevant for magma extraction, we assume H₂O saturation (i.e. H₂O activity of 1) for all calculations by choosing H₂O = 10 wt % (see Ghiorso & Gualda (2015); Gualda & Ghiorso (2014); Harmon et al. (2018) for full discussion). While the choice of using a very large H₂O concentration (i.e. 10 wt % H₂O) renders some simulated magmas fluid oversaturated, H₂O activity is invariably unity, and the presence of exsolved water is immaterial for the calculations. Ghiorso & Gualda (2015); Gualda & Ghiorso (2014) show that the impact of H₂O activity on crystallization pressure is negligible for magmas that are relatively H₂O rich (despite the significant impact of H₂O activity on crystallization temperatures, as summarized in Johannes & Holtz (1996); see also Ghiorso & Gualda (2015)). In other words, the location of the quartz–feldspar cotectic in compositional space is very sensitive to pressure but not to H₂O activity or to temperature. The success of the rhyolite-MELTS geobarometry calculations relies on this experimentally demonstrated behaviour.

We consider phase assemblages including quartz, sanidine, plagioclase, orthopyroxene, and clinopyroxene. As illustrated in Figure 3, only phases that are saturated at the liquidus are in equilibrium with the whole-rock melt composition. The model outputs were manually inspected to determine the liquidus phases. None of our models result in sanidine or clinopyroxene saturation at the liquidus. Sanidine is not present as a phenocryst phase in any of the dome lavas in this study (Supplementary Material B) and it is generally scarce in the central TVZ (Ewart, 1967); it is therefore very unlikely that it is part of the mush mineral assemblage. Clinopyroxene is also very rare as a phenocryst phase in the dome lavas (trace amounts at Deer Hill and Hamurana), so we also assume that they are unlikely to be an important phase in the extracted assemblage; we note, however, that the clinopyroxene calibration in rhyolite-MELTS is potentially problematic for rhyolitic compositions (Brugman & Till, 2019; Seropian et al., 2021). Amphibole and biotite are not modelled appropriately by rhyolite-MELTS (see Gualda et al. (2019b) and our Discussion section for a discussion of implications to this case-study). Although Fe–Ti oxides are included in the rhyolite-MELTS calculations, we do not include them in the final pressure calculations due to their extreme sensitivity to f_O2 (Ghiorso & Evans, 2008). We do not expect Fe–Ti oxides to substantially affect the stability of the other phases we consider here, due to their relatively low abundance and restricted composition spanning only Fe and Ti. We thus consider quartz, plagioclase, and orthopyroxene as possible phases present at the liquidus, which leads to three possible extraction assemblages: quartz + plagioclase, plagioclase + orthopyroxene, or plagioclase + orthopyroxene + quartz (see Pamukçu et al. (2021) for full discussion). Each calculation, therefore, returns two results: the pressure of equilibration between melt and minerals; and the mineral assemblage present at this pressure. For assemblages, including orthopyroxene, the resulting pressure varies with f_O2. It thus results that each composition can yield a range of pressures and assemblages, given the range of f_O2 values we consider.

Monte Carlo simulations

The largest source of random error in the rhyolite-MELTS geobarometer is the analytical error of the melt compositions (Gualda & Ghiorso, 2014; Pamukçu et al., 2021; Pitcher et al., 2021). To estimate the uncertainty of our results, we follow the method of Gualda & Ghiorso (2014) and Pamukçu et al. (2021). We perform Monte Carlo simulations by generating sets of 200 synthetic compositions around known compositions. We selected three representative compositions for analysis (Supplementary Material D): sample Brown_MP13, which gives a plagioclase + quartz pressure of 158 MPa at f_O2 constrained at ΔQFM +1.5; sample Graham_M474, which gives a plagioclase + orthopyroxene + quartz pressure of 292 MPa at f_O2 constrained at QFM; and sample Bowyer_78, which gives a plagioclase + orthopyroxene pressure of 226 MPa at f_O2 constrained at ΔQFM +1. To generate the distribution of synthetic compositions, we use the 1σ error values reported by Weaver et al. (1990) for the University of Canterbury XRF laboratory as referenced by Beresford (1997); Bowyer (2002); Brown (1994); Karhunen (1993); Richnow (1999). We run rhyolite-MELTS simulations on all 600 synthetic compositions following the method described above. From the distribution of results, we calculate statistical parameters (e.g. mean, standard deviation, etc.) for the resulting pressures.

Zircon saturation thermometry

In the same way that we define extraction pressure as the pressure the whole-rock composition equilibrated with the mush mineral assemblage (Figure 1; Gualda et al., 2019b), here we define ‘extraction temperature’ as the temperature the whole-rock magma composition equilibrated with the mush mineral assemblage. We do not use rhyolite-MELTS to calculate the temperature of co-saturation (i.e. the geothermometer equivalent to the rhyolite-MELTS geobarometer) due to the large effect the activity of H₂O can have on liquidus temperatures (Johannes & Holtz, 1996), as well as possible inaccuracies in rhyolite-MELTS temperatures (Gualda et al., 2012; Gardner et al., 2014). Instead, we use zircon saturation geothermometry to assess extraction temperatures. Here, we use the whole-rock Zr concentration to model the temperature the extracted melt equilibrated with a zircon-bearing mush (see Discussion section for discussion of this assumption). We applied both the zircon geothermometer of Watson & Harrison (1983) and the geothermometer of Boehnke et al. (2013) to our dataset of dome compositions.

RESULTS

Whole-rock compositions

Lavas are susceptible to devitrification and alteration due to their slow cooling time. To test the effects of devitrification and alteration on our whole-rock compositions, we consider the compositional variability as a function of the extent of groundmass devitrification reported by the original authors (Figure 4). Inspection of Figure 4 shows that there is no obvious correlation between extent of alteration and composition, and glassy samples come from throughout the study area. These glassy samples span the full compositional range of SiO₂ and K₂O. Some samples with devitrified groundmass have lower Na₂O, which could suggest migration of mobile Na. However, there is no systematic trend as to which samples have lower Na₂O and most of the glassy and devitrified samples overlap. The devitrified samples tend to have lower loss on ignition (LOI) compared to the glassy samples, consistent with dehydration of the groundmass during devitrification. There is no systematic trend between LOI and the alkalis. Overall, these data suggest the dome compositions have been minimally affected by alteration and devitrification.

The compositions of the dome samples can be divided into two groups by their silica content: a high-SiO₂ group with 77.1 to 77.8 wt % SiO₂; and a lower SiO₂ group with 72.4 to 76.5 wt % SiO₂ (Figure 5). The high-SiO₂ group also has lower FeO_T (0.9–1.3 wt %), Al₂O₃ (12.4–12.9 wt %), CaO (0.6–1.2 wt %), and higher K₂O (3.5–4.1 wt %) compared to the low-SiO₂ group (1.2–2.3 wt % FeO_T; 12.7–15.1 wt % Al₂O₃; 0.9–2.1 wt % CaO; 2.7–3.8 wt % K₂O). The high-SiO₂ group samples are predominantly from amphibole-bearing type R1 and R1 + 2 domes (Figure 5). Almost all samples from the NWDC around the Whakamaru caldera have high SiO₂, as well as some domes from the WDC, Rotorua caldera, Reporoa caldera, Ōkataina caldera complex, and the Kapenga volcanic–tectonic depression.

The low-SiO₂ group has a wider range of compositions than the high-SiO₂ group. This group can be subdivided into the more R1-type compositions, which have higher MgO, higher CaO, higher Sr, lower Zr, and lower Y; and the more R2-type compositions (Figures 5 & 6). The R1 domes have plagioclase + quartz + amphibole ± biotite ± orthopyroxene + Fe–Ti oxides phenocryst assemblages, whereas the R2 domes have plagioclase + orthopyroxene ± quartz + Fe–Ti oxides phenocryst assemblages (Supplementary Material B). The R1 + 2 domes have whole-rock compositions that fall between the R1 and R2 groups, and variable phenocryst assemblages of plagioclase ± quartz ± amphibole ± orthopyroxene ± biotite + Fe–Ti oxides. The domes inside the Whakamaru caldera (Maroa dome complex) are type R2, with rare R1 + 2, as are the domes inside the Rotorua caldera. The domes inside the Kapenga volcanic–tectonic depression and around Ōkataina are predominantly R1 and R1 + 2 (Figure 7).

MELTS geobarometry results

Pressure at different QFM

Whilst TVZ rhyolites crystallized under f_O2 conditions ranging from approximately the QFM buffer (ΔQFM +0) to two log units above the QFM buffer (ΔQFM +2) (Deering et al., 2010; Ghiorso & Gualda, 2013), each of the rhyolite magma endmembers has a more restricted range of f_O2 values (R1: ΔQFM +1.5 ± 0.3 1σ; R1 + 2: ΔQFM +0.8 ± 0.3 1σ; R2: ΔQFM +0.5 ± 0.3 1σ; Deering et al., 2010). We therefore refine our results by only considering the calculations at the best estimate of f_O2 for that rhyolite endmember (Figure 8). For R1 domes, we examine ΔQFM +1.0, ΔQFM +1.5, and ΔQFM +2. For R1 + 2 domes, we examine ΔQFM +0.5, ΔQFM +1.0, and ΔQFM +1.5. For R2 domes, we examine ΔQFM +0, ΔQFM +0.5, and ΔQFM +1.0.

The R1 compositions return dominantly quartz + plagioclase solutions at all f_O2, with more plagioclase + orthopyroxene at lower f_O2 (at ΔQFM +1: 49 quartz + plagioclase solutions; 22 plagioclase + orthopyroxene solutions). At all plausible f_O2 values, the R1 domes have very similar distributions: a mode at 150 to 200 MPa and a mode at 250 to 350 MPa (Figure 8a). In contrast, the R2 compositions are deeper, with a single strong mode at 250 to 300 MPa at all f_O2 (Figure 8c). With high f_O2 (ΔQFM +1.0), the R2 compositions return only quartz + plagioclase solutions (n = 96). However, with low f_O2 (ΔQFM +0), R2 compositions return a mixture of quartz + plagioclase (n = 40), plagioclase + orthopyroxene (n = 27), and plagioclase + orthopyroxene + quartz solutions (n = 21). The pressure results for R1 + 2 compositions have the most variability with f_O2 (Figure 8b). At high f_O2 (ΔQFM +1.5), the results are entirely quartz + plagioclase solutions (n = 66); at moderate f_O2 (ΔQFM +1), the modelled mineral assemblage is mostly quartz + plagioclase (n = 58), with some (n = 8) plagioclase + orthopyroxene solutions, and one plagioclase + orthopyroxene + quartz solution. At these moderate to high f_O2, R1 + 2 results yield a wide range of extraction pressures (96–422 MPa) with two weak modes at 150 to 200 MPa and 250 to 300 MPa, similar to the R1 compositions. At low f_O2 (ΔQFM +0.5), R1 + 2 pressures are slightly shallower (97 to 402 MPa), with more plagioclase + orthopyroxene (n = 27) and plagioclase + orthopyroxene + quartz (n = 7) results (quartz + plagioclase n = 36).

The pressure trends for the entire data set are insensitive to the specific choice of f_O2. In the R1 and R2 domes, f_O2 makes very little difference in the distribution of pressure results; but, as expected, f_O2 affects whether orthopyroxene is present in the modelled mineral assemblage (Figure 8). The choice of f_O2 does make a minor difference for the R1 + 2 extraction pressure distribution, at higher f_O2 the pressure distribution is deeper. Individual plagioclase + orthopyroxene pressures are affected by f_O2, with a pressure increase of ~70 MPa per 0.5 log f_O2 units (see example in Figure 3). However, only a relatively small number of samples (n = 27) are affected when considering the total number of compositions in our dataset (n = 239). We conclude that the effect of f_O2 on the overall results is minimal. Given that the results are generally insensitive to f_O2, for the remainder of the discussion, we consider the f_O2 estimated by Deering et al. (2010) as our ‘preferred’ f_O2 value. We then compare results for each magma type at their preferred f_O2. We avoid comparing individual data points, acknowledging that the choice of f_O2 can change the extraction pressure result for a single composition.

Spatial and temporal pressure trends

Across the entire central TVZ region, there are two pressure modes: at 150 to 175 MPa and 250 to 325 MPa, which are equivalent to 5.7 to 6.6 km and 9.4 to 12.3 km depth, respectively (Figure 9). Extraction pressure correlates with magma type: the R2 magmas are predominantly deeper, whilst R1 magmas are shallower; R1 + 2 magmas are evenly distributed between both the shallow and deep modes. Most of the extraction results have quartz in the modelled mineral assemblage. As pressure controls the position of the quartz–feldspar cotectic there is, therefore, a strong correlation between extraction pressure and SiO₂ content (Figure 10) (Tuttle & Bowen, 1958; Johannes & Holtz, 1996; Blundy & Cashman, 2001; Gualda & Ghiorso, 2013b). When our compositions are projected on the haplogranitic ternary of Blundy & Cashman (2001), the projected pressure distribution is qualitatively similar to our rhyolite-MELTS geobarometry results (Figure 10). This illustrates the effect of the shifting quartz–feldspar cotectic with pressure, which rhyolite-MELTS replicates well.

The domes associated with the Whakamaru caldera have two extraction pressure modes strongly correlated with magma type (Figure 9b). The NWDC on the inferred boundary and outside the Whakamaru caldera has shallow (150–200 MPa) extraction of R1 and R1 + 2 magmas (Figure 7). In contrast, the Maroa dome complex within the Whakamaru caldera has deep (250–300 MPa) extraction of R2 magma types, with rare R1 + 2 magmas. Within the limitations of the currently available age dates, none of the dome complexes displays temporal trends in major-element composition, magma type, or extraction pressure (Figure 11). The WDC has the most variability in extraction pressure and composition; the pressures generally get deeper towards the south (Figure 7).

For comparison with the Whakamaru caldera system, Figure 9c shows the pressure results from the domes in the northern and central part of our study area. These also have two pressure modes at 150 to 175 MPa and 250 to 325 MPa; unlike the domes associated with the Whakamaru caldera, the correlation between extraction pressure and magma type is very weak (Figure 9c). The 250- to 325-MPa mode has slightly more R2 domes and the 150- to 175-MPa population has more R1 and R1 + 2 magmas. Generally, the domes within the Ohakuri, Rotorua, and Kapenga caldera structures are extracted from the deeper mode, whilst the domes on the boundaries of the calderas have more varied extraction pressures (Figure 7). An abrupt transition in magma chemistry followed the Whakamaru Group eruptions; the pre-flare-up domes are all R1 and R1 + 2, while R2 magmas only erupt after the Whakamaru Group eruptions (Figure 11).

Monte Carlo analysis

The Monte Carlo simulation for sample Brown_MP13 (quartz + plagioclase assemblage) yields a pressure distribution with mean of 165 MPa and standard deviation of 13 MPa. This estimate is consistent with Pamukçu et al. (2021) and Smithies et al. (2023), who estimate 1σ errors for quartz + plagioclase rhyolite-MELTS geobarometry of 10 to 20 MPa. Our 1σ errors are smaller than Gualda & Ghiorso’s (2014) estimate for the quartz + plagioclase geobarometer of ~25 MPa. Gualda & Ghiorso (2014) estimated errors from a WDS-EMP analysis by Anderson et al. (2000), which had much higher relative errors than the XRF analyses reported here. In particular, Na₂O contributed the most to the uncertainty reported by Gualda & Ghiorso (2014); they use a relative error of 13% for Na₂O, compared to 2% from this study.

For sample Bowyer_78 (plagioclase + orthopyroxene assemblage), 196 of the 200 synthetic compositions return pressure results, with a mean of 228 MPa and standard deviation of 18 MPa. The slightly larger uncertainty associated with the plagioclase + orthopyroxene pressures compared to the quartz + plagioclase case is a result of the effect of f_O2 on orthopyroxene stability (see Pamukçu et al., 2021). Our estimated error for the plagioclase + orthopyroxene geobarometer (1σ of 18 MPa) is smaller than what was estimated by Harmon et al. (2018) (1σ of 26 MPa). Harmon et al. (2018) use glass compositions collected by EDS-SEM, with most of the variability in pressures stemming from the relatively high uncertainty in CaO (relative error of 5%). In contrast, the uncertainty associated with CaO in the XRF data used here is substantially smaller (0.5%), which probably explains our lower resulting uncertainty in pressure values.

For sample Graham_M474 (plagioclase + orthopyroxene + quartz assemblage), 101 of the 200 compositions return a plagioclase + orthopyroxene + quartz result, with a mean pressure of 297 MPa and standard deviation of 12 MPa.

Taken together, these results suggest uncertainties in pressure due to random error are on the order of 15 MPa (1σ), irrespective of the mineral assemblage.

Zircon saturation thermometry results

The extraction temperature results from the calibration of Boehnke et al. (2013) range from 698°C to 830°C, excluding three outlier points at 514°C, 660°C, and 870°C (Figure 12). The results from the Watson & Harrison (1983) calibration are ~50°C warmer, ranging from 750°C to 853°C, excluding three outlier points at 584°C, 717°C, and 884°C (Figure 12). The results from the Watson & Harrison (1983) geothermometer are in better agreement with pre-eruptive temperatures of 750–850°C for other TVZ rhyolites estimated from Fe–Ti oxides and mineral-melt geothermometers (e.g. Schmitz & Smith, 2004; Wilson et al., 2006; Shane et al., 2007; Deering et al., 2010; Allan et al., 2012; Shane & Smith, 2013; Bégué et al., 2014a; Barker et al., 2015; Allan et al., 2017). For the remainder of this discussion, we will therefore focus on the extraction temperature results from the Watson & Harrison (1983) calibration. We emphasize that use of the Boehnke et al. (2013) calibration would only change the absolute temperatures obtained, but not the relative differences, which are more important for us. These zircon-saturation temperatures are independent of the rhyolite-MELTS pressure calculations, which use major-element compositions only.

The temperature distribution is bimodal, with a dominant mode at ~810°C and a secondary mode at ~770°C. The warmer ~810°C mode is predominantly derived from samples from R2 and R1 + 2 domes, whilst the cooler ~770°C mode is entirely resulting from samples from R1 and R1 + 2 domes. Samples from the Whakamaru caldera have strong spatial clustering: samples from the NWDC are predominantly from the cooler mode, whilst samples from the WDC and Maroa are warmer (Figure 13). In other parts of the central TVZ, the extraction temperatures are generally from the warmer mode, and there are no strong spatial trends (Figure 13).

DISCUSSION

What is the depth of magma extraction from the mush?

To infer the depth of the mush from our extraction pressures, we first need to test our assumption that the pressure of equilibration between whole-rock compositions and the modelled mineral assemblages represents extraction. By using whole-rock compositions, we assume perfect extraction whereby the melt does not entrain any crystals from the mush (Figure 1). The bulk rock undoubtedly contains some antecrystic and xenocrystic crystals and crystal cores (e.g. Hildreth, 2004; Wilson et al., 2007; Charlier et al., 2008; Claiborne et al., 2010; Allan et al., 2013; Cooper & Kent, 2014; Sas et al., 2021). However, we argue that the proportion of these components is volumetrically minor and would not shift the bulk composition significantly. For example, antecrystic plagioclase cores from the crystal-poor (~10 wt % crystals) Oruanui eruption make up 0.2% of the bulk volume of the magma (Charlier et al., 2008), which is undetectable in the whole-rock major-element compositions.

Most of the dome lavas in this study are crystal poor (from the domes with reported crystal contents, 89% of the samples are from domes with <20 vol % crystals), thus the melt composition dominates. Within the crystal-poor samples, there is no strong correlation between the reported crystal content of rocks that form each dome and the pressure results for each compositional group (NWDC, other R1 and R2; Figure 14). The samples from domes with >20 vol % are almost all high pressure (300–420 MPa), which could be due to an antecryst population. Crystal-rich lavas containing complex crystal cargo are unlikely to be suitable to calculate extraction pressures. However, these crystal-rich domes are mostly from the same geographic area (north-west Kapenga), so the cluster of high pressures could be a geographic feature. If we exclude the crystal-rich domes, our overall results do not change. The overall range of pressures is from 100 to 375 MPa, and there are two modes at 150 to 200 MPa and 250 to 325 MPa (Figure 14)—in this sense, our results are not changed by the inclusion or exclusion of the more crystal-rich samples, which supports the conclusions presented below.

It is also worth noting that rhyolite-MELTS has been shown to be effective at filtering spurious compositions (Gualda & Ghiorso, 2014; Bégué et al., 2014b; Pamukçu et al., 2021). We expect that a bulk-rock composition with a high proportion of antecrysts or xenocrysts would not return a pressure result as the mixture between melt and arbitrary crystal proportions would not match natural rhyolite compositions suitable for rhyolite-MELTS pressure calculations—for instance, the bulk-rock compositions from the Tauhara dome complex, which is known to include mixing and mingling (Millet et al., 2014), return no pressure results. This suggests that rhyolite-MELTS is effective at weeding out whole-rock compositions that do not represent extracted melts.

The most striking feature of our extraction pressure results is that there are two depths of melt extraction: 150–175 MPa (~6 km) and 250–325 MPa (9–12 km). These modes are quite evident, and they persist even if subsets of the data are excluded (see above). Note that our dataset includes over one hundred eruptions, and the persistence of the two modes suggests that we have a representative dataset of whole-rock compositions from the dome populations. The deeper 9- to 12-km mode correlates with the mid to lower crustal zone (~8–15 km) of partial melt evident in geophysical studies of the modern TVZ (Figure 15). Not only does this agreement give us confidence in our results, but it also suggests that magma mush consistently develops between 9 and 12 km from the present day back to at least 650 ka. It is important to note that the eruptible magma is not necessarily stored in contact with the mush; it may move to a physically separate, shallower, pre-eruptive storage zone (Gualda et al., 2019b; Smithies et al., 2023). Explosive eruptions in the central TVZ are typically stored at ~6 km (~150 MPa), shallower than the deep mush zone revealed here (Bégué et al., 2014b). For the dome eruptions, we cannot determine storage pressures by rhyolite-MELTS geobarometry (see above), so it is currently unknown if dome magmas were stored at shallower levels like many of the caldera-forming eruption magmas, or if they erupted from levels similar to those of extraction.

A smaller number of the extraction pressures are clustered at ~6 km. This overlaps with the ‘typical’ pre-eruptive storage depth of the eruptible magmas feeding most central TVZ eruptions (Figure 15) (Bégué et al., 2014b; Gualda et al., 2018). The shallow extraction pressures may represent mush that is in contact with stored eruptible magma bodies. The dome eruptions may also be tapping residual mush left in the shallow crust following the caldera-forming eruptions. Alternatively, there may be a rheological reason for ascending magmas to stall at this depth forming both magma mush and ephemeral eruptible magma bodies. The ~6-km ‘sweet spot’ for magma stalling is coincident with the brittle–ductile transition zone in the central TVZ (Figure 15; Bryan et al., 1999).

The wide range of extraction pressures raises the question: do individual dome eruptions erupt magmas extracted from a vertically extensive depth range, or are the magmas extracted from a more vertically restricted mush body? It is well established that very large (> > 10 km³) rhyolitic eruptions can contain multiple magma types which could be extracted from different magma mush bodies, stored in separate magma bodies, and only mingled during eruption (Cooper et al., 2012; Gualda & Ghiorso, 2013a; Cashman & Giordano, 2014; Bégué et al., 2014a; Myers et al., 2016; Kennedy et al., 2018; Swallow et al., 2018; Swallow et al., 2019; Pearce et al., 2020; Gualda et al., 2022). In the central TVZ, there is evidence that smaller (<10 km³) eruption episodes can contain multiple magma batches with contrasting crystallization pressures, such as the <21.8 ka Tarawera volcanic complex in the Ōkataina volcanic centre (Shane et al., 2007; Shane et al., 2008a; Shane et al., 2008b; Storm et al., 2012). To investigate the potential for multiple magma types within single domes, we examine in greater detail the 10 domes for which we have more than five compositions (Figure 16). Three domes (Maungawhakamana, Te Tarata, Wereta) have variable major- and trace-element compositions (Figure 16a). Both Maungawhakamana and Te Tarata have high variance in the extraction pressures (Figure 16b). This suggests that these domes are comprised of multiple magmas extracted from mush at different depths, erupted either during the same event (i.e. multiple magma batches) or during different eruptive episodes. The remaining eight domes have very narrow ranges in extraction pressures, with standard deviations of 14 to 33 MPa. These small variances are similar to the <20 MPa 1σ errors we estimate from analytical uncertainty. This gives us confidence that we can assign extraction pressures to a dome eruption with a high degree of certainty. It also suggests that—in many cases—the domes are one magma batch extracted from a single mush body with limited vertical extent.

What is the mineral assemblage in the mush?

A major uncertainty in our extraction models is that we do not know a priori what the mineral assemblage of the mush is, so we do not know with what assemblage our extracted melt equilibrated. We follow the approach of Gualda et al. (2019b), who consider two possible mush mineral assemblages in the central TVZ: a quartz-bearing ‘granitoid’ assemblage with quartz + plagioclase ± orthopyroxene ± amphibole ± biotite (quartz + plagioclase and plagioclase + orthopyroxene + quartz results); or a quartz-absent ‘diorite’ assemblage with plagioclase + orthopyroxene ± amphibole ± biotite (plagioclase + orthopyroxene results). There is evidence for both assemblages from the mineralogy of erupted plutonic lithics in the central TVZ (Browne et al., 1992; Burt et al., 1998; Brown et al., 1998a; Deering et al., 2011b; Graeter et al., 2015).

Rhyolite-MELTS does not have good models for amphibole or biotite, due to the complexity of these phases and the lack of experimental data to calibrate them. It is likely that amphibole was part of the mush mineral assemblage in many of the R1 magmas, particularly in the Ōkataina caldera (Burt et al., 1998; Deering et al., 2011b). Importantly, these are the magmas most likely extracted from quartz-bearing sources, thus the quartz + plagioclase geobarometer is still applicable. The crystallization of amphibole is unlikely to affect the calculated saturation pressure or temperature of quartz, as partitioning of SiO₂ into quartz is much higher than into amphibole, and amphibole abundances are small even in the most crystal-rich magmas. Similarly, plagioclase saturation will not be significantly affected by amphibole as the partitioning of CaO, Na₂O, and SiO₂ into plagioclase is much higher than into amphibole. We therefore expect our quartz + plagioclase pressures to be unaffected by the presence of amphibole. Amphibole may affect orthopyroxene saturation more significantly, depending on the composition of the amphibole and the orthopyroxene. Most of our pressure results are quartz + plagioclase, so this effect is not significant.

In the absence of external constraints, the rhyolite-MELTS geobarometer can be used to evaluate the suitability of different extraction assemblages (Figure 3; Gualda et al., 2019b; Pamukçu et al., 2021). It is an interesting feature of our results that the quartz-bearing results dominate, suggesting that most rhyolites are extracted from a quartz-bearing mush (Figure 8). Quartz-bearing plutonic lithics are very common in the central TVZ (Burt et al., 1998; Brown et al., 1998a; Graeter et al., 2015). There is strong evidence that quartz crystallizes down to at least the mid crust, as many TVZ rhyolites have quartz-hosted melt inclusions with maximum H₂O-CO₂ pressures greater than 200 MPa (Shane et al., 2007; Shane et al., 2008a; Smith et al., 2010; Johnson et al., 2011). The presence of quartz also offers an explanation to some of the unusual features of the whole-rock compositions. For example, the SiO₂ content of the NWDC is very tightly constrained (mean SiO₂ 77.3 wt %, 1σ = 0.5) despite the domes having erupted over a period of at least 60 kyrs. In rhyolites, silica content is controlled by cotectic crystallization of quartz and plagioclase (Tuttle & Bowen, 1958; Johannes & Holtz, 1996; Blundy & Cashman, 2001; Gualda & Ghiorso, 2013b), so the tightly constrained SiO₂ content reflects crystallization at a very narrow of depths. We can therefore offer a simple, physical explanation for an otherwise surprising geochemical result: magmas that erupted to form the NWDC domes were repeatedly extracted from a mush body at ~6 km depth.

There is a distinct clustering of SiO₂ content between 74–76.5 wt % and 77–77.5 wt % and correlated clusters in FeO_T and K₂O (Figure 5). In a quartz–feldspar system, the SiO₂ content is strongly controlled by the pressure of crystallization and the position of the quartz–feldspar cotectic (Figure 10; Blundy & Cashman, 2001; Gualda & Ghiorso, 2013b; Johannes & Holtz, 1996; Tuttle & Bowen, 1958), thus the two modes in our extraction pressure are correlated with the two modes in SiO₂ contents (Figure 10). The clusters in other elements (e.g. FeO_T) suggest that other phases (e.g. orthopyroxene) also present and equilibrate with the melt in the mush, showing that the mush is multiply saturated (see Blundy, 2022). Our model offers a mechanism to explain the compositional trends, with the ‘gap’ in SiO₂ reflecting contrasting magma mush depths.

What is the temperature of magma extraction from the mush?

Here we use zircon saturation geothermometry to examine relative differences in extraction temperature. The choice of calibration for the zircon saturation geothermometer changes the absolute temperature values, but both calibrations show clear differences in relative temperature between the domes (Figure 12). The zircon saturation temperatures are clustered in two groups, a cooler mode at ~770°C and a warmer mode at ~810°C. These temperature differences are due to absolute variation in the Zr content of the magmas. We consider three explanations for the distribution of Zr and therefore temperature: 1) zircon undersaturation, which would make the temperatures minima; 2) zircon inheritance, which would make the temperatures maxima; 3) true temperature variation between mush bodies.

For the dome magmas, zircon undersaturation is unlikely. Accessory zircons are reported by Bowyer (2002) in rhyolite lavas of the Ōkataina, Rotorua, and Kapenga volcanic centres. Many of the original sources of the data in this study do not specify zircon presence or absence (Ewart, 1967; Brown, 1994; Leonard, 2003). Zircon is present in many other TVZ rhyolites (Brown & Fletcher, 1999; Charlier et al., 2003; Charlier et al., 2005; Shane et al., 2005; Charlier et al., 2010; Barker et al., 2014; Cooper et al., 2014; Rubin et al., 2016). In the silicic rhyolites (>70 wt % SiO₂) in this study zircon is very likely to be saturated.

The antithetical problem to zircon undersaturation is zircon inheritance, whereby additional zircon crystals are entrained by the extracted magma. The effect is to raise the Zr concentration in the whole-rock composition, thus artificially increasing the zircon saturation temperatures (Miller et al. (2003), see also correlation between Zr and temperature in Figure 12). Zircon inherited from the crust and from previous magmatic episodes are common in rhyolitic eruptions from the TVZ, although the inherited zircons are volumetrically less prevalent than phenocrystic zircons (Charlier et al., 2003; Charlier et al., 2005; Charlier et al., 2010; Rubin et al., 2017). Figure 12 compares the crystal content of the domes to the zircon saturation temperatures. The low crystal content domes should have relatively little zircon inheritance. Notably, in the very low crystal content (<5%) lavas, there is a bimodal spread in temperature from 750°C to 814°C. Variable temperatures between lavas of the same crystal content implies that the zircon saturation temperature of these lavas is not controlled by zircon inheritance, but by the true temperature of crystallization in the mush.

Within the low crystal content lavas from the Whakamaru caldera there is a positive correlation between extraction pressure and temperature (Figure 12). The R1 domes on the caldera margins (NWDC) have both lower extraction pressure and lower extraction temperature than the intra-caldera R2 domes (Maroa) (Figure 13). The solidus curve of a water-saturated haplogranitic system predicts the opposite pressure–temperature relationship (see Fig. 2.13 of Johannes & Holtz (1996)). To account for this, the magmas must have different water activities. If the R2 magmas have a lower water activity than the R1 magmas, their temperature would be higher at higher pressure. This pressure–temperature relationship therefore confirms that the deeper, hotter, R2, orthopyroxene-bearing post-caldera Maroa rhyolites in the Whakamaru caldera are drier than the shallower, colder, amphibole-bearing caldera-bounding NWDC rhyolites.

Adjustment of the caldera magmatic system following a ‘super-eruption’: Whakamaru case study

The Whakamaru caldera is an excellent case-study for the state of the magmatic system following a ‘super-eruption’-scale event in the TVZ. We focus on the abundant post-caldera domes to investigate the rearrangement of the magmatic system following caldera formation.

The extraction pressures and temperatures for domes associated with the Whakamaru caldera show distinct spatial clustering, suggesting they tapped two discrete regions of the magmatic system. The voluminous intra-caldera Maroa dome magmas were extracted from ~225 to ~350 MPa, overlapping with the extraction depth of the Whakamaru ignimbrite magmas (Figure 11a; Smithies et al., 2023). However, the Maroa domes are characterized by R2 compositions, contrasting with the R1 Whakamaru Group compositions (Figure 11a; Deering et al., 2010). Eruption of hotter, drier, more reducing magma after caldera-forming eruptions is a common feature of rhyolitic magma systems (e.g. Kos-Nisyros volcanic centre, Bachmann et al., 2012; Valles caldera, Wilcock et al., 2013; Long Valley caldera, Hildreth et al., 2017). This compositional shift is often attributed to new mafic melts recharging the magmatic system (Bouvet de Maisonneuve et al., 2021). If this is the case in Maroa, recharge must have been voluminous and must have rapidly evolved to generate ~600 km³ of rhyolitic Maroa magmas within 45 kyrs of the Whakamaru Group eruptions (Leonard, 2003).

In contrast to Maroa, the NWDC on the north-western boundary of the Whakamaru caldera has similar R1 compositions to the Whakamaru ignimbrite (Figure 11). The magmas that compose the NWDC were extracted from shallow depths of 4 to 7 km (100–200 MPa). Although they have typical R1 amphibole-bearing mineral assemblages, they are unusually crystal-poor for R1 magmas (<8 vol. %, cf. 10–45 vol. % for most R1; Deering et al., 2010), they have low extraction temperature (~770°C), and they are also some of the most evolved compositions within our study area (Figure 5). The depth of this extra-caldera reservoir overlaps with pre-eruptive storage depths of the Whakamaru ignimbrite magmas (Figure 11). This, in addition to the highly evolved compositions, may indicate that the mush body(s) feeding these domes are the partially crystalline remnants of the pre-eruptive storage region that fed the Whakamaru ignimbrites. Unpublished ⁸⁷Sr/⁸⁶Sr ratios from the NWDC are similar to the Whakamaru ignimbrite magmas, suggesting a similar source (Brown, 1994). Further work is needed to establish the relationship between the NWDC and the older Whakamaru ignimbrite reservoir; for example, zircon age spectra could reveal whether the NWDC mush is entirely new magma (e.g. the post-Oruanui eruptions; Barker et al., 2014), or whether it includes zircon (and melt) from the old Whakamaru system. In any case, it is interesting to note that remains of the Whakamaru system may have persisted after the main group of eruptions, but that they also may have coexisted laterally with the R2 magmatic systems that fed the Maroa eruptions within the caldera structure.

The contrast between the Maroa and the NWDC domes highlights the spatial heterogeneity of the magmatic system from which they derive (Figure 17a). Interestingly, it mirrors the spatial heterogeneity of the magmatic system supplying the Whakamaru caldera-forming eruptions, which were fed by multiple, laterally adjacent, compositionally distinct magma bodies (Brown et al., 1998b; Harmon, 2022; Smithies et al., 2023; Harmon et al., 2024). The intra-caldera Maroa dome complex is strongly influenced by extensional tectonics (Wilson et al., 1986); the domes are aligned along the dominant tectonic lineament in the main rift axis and several domes were fed by magmatic dyking aligned along faults (Leonard, 2003). Rifting may have enabled extraction of the voluminous dome magmas from the mid to lower crust without stalling in a shallower mush zone. Unlike many calderas globally (e.g. Lake City, Hon et al., 1989; La Pacana, Gardeweg & Ramírez, 1987; Long Valley, Hildreth et al., 2017; Toba, de Silva et al., 2015; Mucek et al., 2017; Mucek et al., 2021; Valles, Kennedy et al., 2012; Wilcock et al., 2013), the Whakamaru system shows no evidence for post-collapse resurgence of the caldera floor caused by new magma intruding below the caldera. Voluminous eruption in Maroa appears to have accommodated the recharging magmas without requiring structural uplift. In contrast to Maroa, the NWDC and WDB are outside of the actively rifting region, and they are generally aligned along the caldera margin mapped by gravity surveys and the associated ring faults (Wilson et al., 1986; Stagpoole et al., 2020). Post-caldera domes commonly use the pathways created by the caldera-collapse structures (e.g. Campi Flegrei, Forni et al., 2018; La Pacana, Gardeweg & Ramírez, 1987; Long Valley, Hildreth, 2004; Valles, Wilcock et al., 2013). The caldera boundary zone may have provided a pathway for the cool, shallow, NWDC magmas to erupt rather than fossilize as intrusions (Loewen et al., 2017).

The regional magmatic system

The caldera-forming eruptions of the ignimbrite flare-up erupted ~3000 km³ of magma (DRE). The dome eruptions represent at least another 740 km³ of magma. The dome eruptions and the magma systems feeding them has previously been under-represented (cf.Gualda et al., 2018 ; Smithies et al., 2023), but together with the caldera-forming eruptions they reveal the regional magma system of the central TVZ before, during, and after the ignimbrite flare-up.

Pre-flare-up

Prior to the start of the ignimbrite flare-up at 350 ka, domes erupted around the northern margins of Ōkataina and the margins of Rotorua and Reporoa calderas (Figure 2b). There is no record of pre-flare-up domes in the central region; these were probably buried and destroyed by caldera collapse, younger faulting, and sedimentation. The pre-flare-up domes are as old as 531 ± 10 ka (Whakapoungakau, Ōkataina), extending through to the beginning of the flare-up (330 ± 1 ka, Matawhaura, Ōkataina) (Leonard et al., 2010). There is also a probable caldera collapse in Ōkataina before this period, the Utu ignimbrite at 557 ± 3 ka (Leonard et al., 2010; Deering et al., 2011b). This volcanism in the northern study area overlaps in time with the range of zircon ages from the Whakamaru ignimbrite, which extend back to 608 ± 20 ka (Brown & Fletcher, 1999), indicating that magmatism occurred throughout the entire central TVZ region in the build-up to the flare-up.

These pre-flare-up eruptions all have R1 and R1 + 2 compositions, similar to the Whakamaru Group ignimbrites that followed them (Figure 11; Deering et al., 2010). These domes are compositionally dissimilar to the caldera-forming eruptions that later erupted through the same calderas, for example, the domes around Rotorua caldera have hornblende-bearing R1 compositions, in contrast with the orthopyroxene-bearing R2 Mamaku ignimbrite magmas erupted from the Rotorua caldera (Milner et al., 2003; Deering et al., 2010; Ashwell et al., 2013). The processes driving R1 magmatism occurred across the entire central TVZ region and were not confined to just the southern Whakamaru region. These pre-flare-up domes may be better treated as precursors to the Whakamaru Group eruptions, as part of the build-up to the flare-up eruptions.

The pre-flare-up extraction pressures have a wide distribution from 125 to 380 MPa. There is no spatial clustering; for example, the domes around Ōkataina have disparate pressure results despite being <10 km apart (Figure 7). This suggests that the dome eruptions were fed by discrete, vertically separated mush bodies within a very small area (Figure 17b). Variable extraction depths may be characteristic of the early phase of caldera magma system growth, before the large caldera-forming mush bodies begin to grow, homogenize, and trap ascending magmas.

Syn-flare-up

The Whakamaru Group eruptions mark an abrupt transition in the compositions and eruptive style of the central TVZ. After the Whakamaru Group eruptions, the entire central TVZ erupts R1 + 2 and R2 magmas (Deering et al., 2010). Both caldera-forming eruptions and domes are dominated by R1 + 2 and R2 compositions (Figure 11), suggesting a shift in either the composition of the mantle source or of the mafic fractionation processes (Deering et al., 2010).

Following the Whakamaru caldera-forming eruption, the subsequent caldera-forming eruptions have magmatic systems that are initially deep and vertically condensed, before become progressively vertically extensive, recording vertical growth and maturation of the regional magma system (Gualda et al., 2018; Smithies et al., 2023) (Figure 11). At the same time as the caldera-forming eruptions tapped this maturing regional magma system, the dome eruptions tapped mush bodies that span the entire 4- to 15-km range of the mid to lower crust (Figure 11). Without better age constraints on the dome eruptions, we are uncertain when they erupted between the caldera-forming eruptions. Many of the domes have extraction pressures that overlap with the extraction pressures of the caldera-forming eruptions (200–400 MPa), indicating that they could have the same source mushes (Figure 15). Individual caldera-forming eruptions tap a wider range of extraction pressure compared to individual dome eruptions, suggesting that the larger eruptions tap larger mush bodies. The caldera-forming extraction pressures predominantly overlap with the deeper 9- to 12-km extraction pressure mode seen in the dome magma systems, rather than the shallower extraction mode at ~6 km seen in the dome magma systems (Figure 15), possibly because the deeper mush is better suited to generating very large volumes of eruptible magma, whereas the shallow crust is mostly occupied by ephemeral storage and smaller, more rapidly solidified mush bodies.

Post-flare-up

Domes immediately following the flare-up (240–150 ka) are concentrated within the Rotorua, Ohakuri, and Kapenga calderas (Figure 2b). The volume of lava from Rotorua and Ohakuri is significantly smaller than Maroa, suggesting they are in the waning stages of the flare-up magma supply. There have been no further caldera-forming eruptions from Whakamaru, Rotorua, Ohakuri, and Reporoa calderas. These domes represent the end of flare-up magmatism and the shifting of volcanism in the modern TVZ towards Ōkataina in the north and Taupō in the south (Barker et al., 2020).

The domes in the northern area of the Kapenga volcanic–tectonic depression contrast with the other post-flare-up domes. Unusually, they have R1 compositions, despite being extracted from 10 to 15 km (Figure 7). The northern Kapenga domes have compositions that contrast with the preceding R2 caldera-forming eruptions sourced from the Kapenga region (Chimp, Pokai, Mamaku, and Ohakuri ignimbrites) (Figure 17c). This group of ‘Kapenga’ domes may instead be related to present R1 and R1 + 2 magmatism in the Ōkataina system rather than the older flare-up system (Deering et al., 2010; Sas et al., 2021).

Summary

Within the central TVZ there are two depth intervals that are repeatedly occupied by magma at ~6 km and ~ 9 to 12 km (Figure 15). These depths are used by both small dome-forming and large caldera-forming eruptions. This raises a mechanistic question—why are these levels of the crust preferential for magma stalling and mush development? Consistent shallow storage depths have been observed in many rhyolitic systems globally (Huber et al., 2019), but with the extraction geobarometer, we can also reveal a deeper level of consistent magma accumulation. The long record of dome eruptions (Figure 11) shows that melt is often (if not always) present beneath the calderas, even if the system is not primed for a caldera-forming eruption. Although the caldera-forming eruptions and the dome eruptions tap mush at the same levels in the crust, more work is needed to confirm the connection—are these different styles of eruptions tapping the same magma bodies? Or are the depths the same because of external controls?

The magma depth that is tapped during the small eruptions is controlled by volcanic and tectonic features (caldera-boundary zones and normal faulting). There is also a correlation between the composition of the domes, the depth of magma extraction, and the location of the domes relative to volcanic features. The small eruptions help us untangle the spatial complexity of the caldera magma systems to a level of detail that is not possible with destructive large caldera-forming eruptions.

CONCLUSIONS

1. The rhyolite-MELTS geobarometer can be used to constrain extraction pressures. Pressure uncertainties due to analytical error are <20 MPa (1σ). The extraction pressure results are correlated with SiO₂ content. Our results suggest central TVZ rhyolites are predominantly extracted from quartz-bearing mush, thus the pressure of quartz–feldspar crystallization at the cotectic controls SiO₂ content. At certain f_O2 values, quartz-free extraction assemblages are possible, but the extraction pressure trends are similar. The geobarometer is readily applied to large datasets, such as our 103 dome eruptions, enabling us to assess regional trends in extraction depths.

2. There are two main magma extraction (mush) ‘zones’ in the central TVZ: one ~6 km (150–175 MPa) and another between 9 and 12 km (250–325 MPa). These zones are persistent throughout the build-up to, during, and after the ignimbrite flare-up, and they are persistent across the central TVZ. The deep zone correlates with previous geophysical and monitoring studies, which show partial melt below ~8 km. The lower to mid crust is a therefore a significant region for the fractionation of rhyolitic magmas that is persistent through time. Active rifting may promote extraction from these deeper mush bodies. The shallower mush zone at ~6 km is subordinate and more transient, correlating with shallow pre-eruptive storage of the caldera-forming magmas. Stalling at these shallow depths is coincident with the brittle–ductile transition.

3. Post caldera-collapse domes in the Whakamaru caldera are extracted from two spatially distributed reservoirs. Within the caldera, Maroa has hotter (~810°C), more reducing, orthopyroxene-bearing magmas extracted from deeper levels (~9–13 km, 225–350 MPa), and represents onset of a new system enhanced by mafic recharge and rifting. On the boundaries of the caldera, the NWDC has very evolved magmas from cooler (~770°C), shallow mush (~4–7 km, 100–200 MPa), and reveals progressive crystallization and ‘death’ of the old system. Both magmatic subsystems erupted at generally the same times, highlighting the influence of structural features, such as tectonic rifting and caldera-bounding faults on magma extraction.

4. Before the ignimbrite flare-up, magmas were generally of the colder, wetter, more oxidizing R1 rhyolite type and were extracted from vertically dispersed mush bodies, hinting at the maturation of caldera-forming magma systems from dispersed magma systems to larger (>10 km³) bodies. During and after the flare-up, extraction and composition trends outside the Whakamaru caldera are complex and are likely due to localized structural and magmatic features.

DATA AVAILABILITY STATEMENT

The data underlying this article are available in the article and in its online supplementary material.

Acknowledgements

We thank Lydia Harmon, Graham Leonard, Chad Deering, Kari Cooper, Adam Kent, and Alex Nichols for helpful discussions. Luca Caricchi and Marco Brenna gave thoughtful comments on an earlier version of this manuscript. We thank Madeleine Humphreys and Georg Zellmer for editorial handling, and Scott Bryan and two anonymous reviewers for their detailed comments. We also acknowledge the many students whose data are cited in this study, and Jim Cole, who supervised most of them.

REFERENCES