# Hydrothermal Cooling as a Requirement for Short Storage of Silicic Magmas

**Authors:** Liam J. Kelly, Guilherme A. R. Gualda, Darren M. Gravley, David E. Dempsey  
**Published in:** Geochemistry, Geophysics, Geosystems (2021)  
**DOI (version of record):** https://doi.org/10.1029/2021gc009794  
**This copy:** Published version. Text extracted from the PDF at https://ddempsey.github.io/papers/2021-kelly-hydrothermal-cooling-requirement-short-storage/kelly2021-hydrothermal-cooling-requirement.pdf; equations, tables and figure placement may be imperfect, so cite the DOI.  
**Notice:** Published version, posted under AGU's author self-archiving policy. © 2021 American Geophysical Union. All Rights Reserved. Version of record: https://doi.org/10.1029/2021gc009794  

---


# **Hydrothermal Cooling as a Requirement for Short Storage of Silicic Magmas** 

## **RESEARCH ARTICLE** 

### 10.1029/2021GC009794 

##### **Key Points:** 

### **Liam J. Kelly**<sup>**1**</sup> **, Guilherme A. R. Gualda**<sup>**1**</sup> **, Darren M. Gravley**<sup>**2**</sup> **, and David E. Dempsey**<sup>**3,4**</sup> 

- In order to cool silicic magma bodies over decades to centuries, hydrothermal systems must drive the heat loss 

1Department of Earth and Environmental Sciences, Vanderbilt University, Nashville, TN, USA, 2School of Earth and Environment, University of Canterbury, Christchurch, New Zealand,<sup>3</sup> Department of Civil and Natural Resources Engineering, University of Canterbury, Christchurch, New Zealand,<sup>4</sup> Department of Engineering Science, University of Auckland, Auckland, New Zealand 

- The current heat output of the Rotorua is not sufficient for cooling a 150 km<sup>3</sup> magma body over decadal timescales 

**Abstract** Large (>50 km<sup>3</sup> ) magma bodies that contribute to caldera-forming eruptions in the Taupō Volcanic Zone (TVZ), New Zealand have shown evidence for crystallization and storage on the order of decades to centuries prior to eruption. We compare various conductive and convective heat extraction models to determine which is more likely to enable the heat loss required for short magma storage. We conclude that purely conductive models only account for heat loss on millennial or longer timescales. We also show that convective hydrothermal systems with heat output of 1,000 MW in magnitude are required for decadal heat extraction from a large, contiguous magma body; yet, heat output similar to present-day conditions would be suitable for cooling magma distributed as a patchwork of smaller magma bodies. This study shows the potential connection with heat flow at the surface and the presence of a magma body, providing future directions for monitoring restless calderas in environments with a wealth of hydrothermal activity like the TVZ. 

- Cooling, silicic magma bodies in the Taupō Volcanic Zone are likely to be influenced by hydrothermal systems, since there are currently ∼20 active systems 

##### **Supporting Information:** 

Supporting Information may be found in the online version of this article. 

##### **Correspondence to:** 

L. J. Kelly, lj.kelly@vanderbilt.edu 

**Plain Language Summary** Previous research has shown that eruptible magma bodies in the Taupō Volcanic Zone (TVZ), New Zealand may only have been stored in the crust for less than 100 years. These magma bodies have contributed to highly explosive, caldera-forming eruptions, so understanding the timescales over which heat is lost from these magma bodies prior to eruption could be beneficial for monitoring these caldera environments. In order to understand how these magma bodies can lose heat this quickly, we model heat loss through conductive and convective methods, taking advantage of the existence of hydrothermal systems that currently reside in the TVZ. We find that hydrothermal systems are required for rapid heat loss over decades to centuries. Knowing the scale of the power output from such a hydrothermal system could be important for understanding how quickly possible eruptible magma can be cooled, and could be combined with measurements determining the presence or absence of magma to assess potential hazards in the TVZ and beyond. 

##### **Citation:** 

Kelly, L. J., Gualda, G. A. R., Gravley, D. M., & Dempsey, D. E. (2021). Hydrothermal cooling as a requirement for short storage of silicic magmas. _Geochemistry, Geophysics, Geosystems_ , _22_ , e2021GC009794. https://doi. org/10.1029/2021GC009794 

Received 15 MAR 2021 Accepted 17 AUG 2021 

##### **Author Contributions:** 

**Conceptualization:** Guilherme A. R. Gualda, Darren M. Gravley **Data curation:** Liam J. Kelly, Guilherme A. R. Gualda **Formal analysis:** Liam J. Kelly **Funding acquisition:** Liam J. Kelly **Investigation:** Liam J. Kelly **Methodology:** Liam J. Kelly, David E. Dempsey **Project Administration:** Darren M. Gravley 

**Resources:** Guilherme A. R. Gualda, Darren M. Gravley **Software:** Liam J. Kelly **Supervision:** Guilherme A. R. Gualda, Darren M. Gravley, David E. Dempsey **Validation:** Liam J. Kelly, Guilherme A. R. Gualda **Visualization:** Liam J. Kelly 

© 2021. American Geophysical Union. All Rights Reserved. 

## **1. Introduction** 

Ignimbrite-producing, caldera-forming eruptions are some of the most intense, explosive eruptions on the planet, occurring in a wide variety of tectonic environments (Bryan et al., 2010; de Silva & Gosnold, 2007; Elston, 1984; Gravley et al., 2016; Lipman, 1984). Understanding the heat flux associated with large caldera-forming eruptions is key since, if there is a significant increase in heat flux of a hydrothermal system that is already in existence when a large silicic magma body is being stored, then the measurement of heat flux at the surface of calderas in a state of unrest may offer insight into the eruptibility of magma below the surface. 

Crystal growth times inferred from geochronometry and geochronology suggest that magma bodies are typically not stored for longer than centuries to a few millennia (Allan et al., 2013; Charlier et al., 2008; K. M. Cooper & Kent, 2014; Druitt et al., 2012; Gualda & Sutton, 2016; Gualda et al., 2018; Gualda, Pamukcu, et al., 2012; Pamukcu, Gualda, Bégué, & Gravley, 2015; Shamloo & Till, 2019). In the Taupō Volcanic Zone (TVZ), New Zealand, crystallization may be even shorter, spanning only decades to a few centuries, as demonstrated by Ti-in-quartz geochronometry (Gualda et al., 2018). This suggests that eruptible magma bodies in this region are only stored for decades to centuries, at relatively constant pressures, prior to catastrophic, caldera-forming eruptions (Gualda et al., 2018; Pamukcu, Gualda, Bégué, & Gravley, 2015). While Ti diffusivities in quartz have been recently questioned (Jollands et al., 2020a, 2020b; see also Audétat 

1 of 17 

KELLY ET AL. 

**Geochemistry, Geophysics, Geosystems** 

10.1029/2021GC009794 

**Writing – original draft:** Liam J. Kelly **Writing – review & editing:** Liam J. Kelly, Guilherme A. R. Gualda, Darren M. Gravley, David E. Dempsey 

et al., 2021), other sources of evidence suggest that these slow diffusivities do not prevail under magmatic conditions (Gualda & Pamukçu, 2020; Pamukcu, Gualda, Bégué, & Gravley, 2015). Here we emphasize that crystallization timescales constrain the magma cooling timescales in the crust, as well as storage timescales of an eruptible magma. 

The question then arises: how is heat extracted so that crystallizing, eruptible magma is only stored for decades to centuries? In order for a magma body to crystallize, it must lose heat, which includes sensible heat loss, which leads to decreasing temperatures, and latent heat loss, which is necessary for crystallization of different mineral phases and bubble exsolution (Perrot, 1998). The complexity of the problem increases if convection within the magma chamber is taken into account; however, a convecting rhyolite can accelerate the cooling of a reservoir by 15% (Andersen & Weis, 2020). While we do not consider the effect of convection in the magma body, any effect would only decrease heat extraction timescales calculated in this study. Heat extraction from magma bodies as a result of conductive heat flow has been previously studied (Annen et al., 2015; Gualda, Pamukcu, et al., 2012), as well as the impact of permeability of surrounding rock on the mode of heat flow (Hayba & Ingebritsen, 1997), the longevity of hydrothermal systems sustained by ultramafic intrusions (Cathles et al., 1997), and characterization of the exchange between hydrothermal fluids and magmatic intrusions (Norton & Taylor, 1979). Convective heat extraction as it relates to geothermal wells, economic viability, temperatures near the liquidus, and different magma types has been investigated, but not the timescales of said heat extraction (Hardee, 1981; Hardee & Dunn, 1981; Hardee & Larson, 1977). 

The mechanisms for heat extraction from a magma body that could contribute to crystallization of these large (>50 km<sup>3</sup> ) magma bodies on the order of decades to centuries are not well understood. In this study, we compare the effects of conductive heat transfer—through contact of the magma body with the overlying crust—and convective heat transfer—caused by the movement of heat by the flow of fluids in the crust overlying the magma body (Boden, 2016). We follow Gualda, Pamukcu, et al. (2012), Pamukcu, Gualda, Ghiorso, et al. (2015), and Gualda et al. (2018) in using rhyolite-MELTS to model the change in enthalpy (total heat loss) during crystallization of a magma body at constant pressure. Then, we compare heat transfer within the overlying crust by various conductive and convective heat flow models to calculate the time required to lose the total heat obtained from rhyolite-MELTS. To test the sensitivity of our results to variability inherent in the model parameters, we used a Monte Carlo approach. 

We constrain our model by using geological information from the Mamaku ignimbrite, the product of the ca. 240 ka Rotorua caldera-forming eruption in the central TVZ (Bégué et al., 2014; Gravley et al., 2007). This region is an ideal place for this research, since it is home to hydrothermal systems sitting within and on the margins of large silicic calderas (Bibby et al., 1995; Bégué et al., 2014; Wilson & Rowland, 2016). There is also a wealth of knowledge regarding geothermal systems due to geothermal energy exploration and production in the region (Bibby et al., 1995; Hedenquist, 1990; Rowland & Simmons, 2012), which has provided us with an understanding of the cross-sectional fluid and rock properties, as well as the reservoir conditions of these dynamic systems. By characterizing the influence of hydrothermal systems on magma bodies and providing a link between crustal processes and magmatic processes, we expand the toolkit for those who monitor and assess the hazards associated with restless calderas. 

### **1.1. Geological Background** 

The central TVZ includes deposits of many caldera-forming eruptions, having erupted ∼6,000 km<sup>3</sup> of rhyolitic material over a period of ∼1.6 Ma (Houghton et al., 1995; Wilson et al., 1995). The TVZ flare up (Gravley et al., 2016) was a time of highly productive rhyolitic volcanism that lasted from ∼350 to 240 ka, with three main pulses that led to a total of ∼3,000 km<sup>3</sup> of tephra erupted from multiple caldera-forming eruptions. 

The Mamaku ignimbrite erupted in the third pulse of the TVZ flare up from the Rotorua Caldera (Figure 1) at ca. 240 ka (Gravley et al., 2007); it has a volume of >145 km<sup>3</sup> and was sourced from multiple magma batches extracted from the same crystal mush (Bégué et al., 2014; Gravley et al., 2016). The Mamaku ignimbrite was chosen for this study because the Rotorua caldera currently hosts a geothermal system, and Ti-in-quartz geochronometry of pumice-hosted quartz in the Mamaku ignimbrite shows crystallization on the order of decades to centuries (Gualda et al., 2018). Despite recent disputes about Ti-in-quartz timescales (Gualda & Pamukçu, 2020; Jollands et al., 2020a, 2020b), crystallization timescales are confirmed by other 

2 of 17 

KELLY ET AL. 

**Geochemistry, Geophysics, Geosystems** 

10.1029/2021GC009794 

<!-- Start of picture text -->
N<br>200 km<br>38 S<br>0 10 20 km<br>38.5 S<br>176 E 176.5 E<br><!-- End of picture text -->

methods (e.g., melt inclusion faceting, see Gualda, Ghiorso, et al., 2012; Pamukcu, Gualda, Bégué, & Gravley, 2015), as well as other independent estimates (see Pamukcu, Gualda, Bégué, & Gravley, 2015). 

Mamaku pumice is high-silica rhyolite with a mineralogy of primarily plagioclase, quartz, and orthopyroxene, and crystal contents of pumice are <10 wt% crystals (Bégué et al., 2014). Storage pressures for the Mamaku ignimbrite range between 50 and 160 MPa (∼1.8–5.5 km depth) (Bégué et al., 2014; Gualda et al., 2018). 

The present-day heat transfer of the Rotorua geothermal system is estimated to be 470 MW over an area of 30 km<sup>2</sup> (Bibby et al., 1995). The broader TVZ is a very active geothermal area, with over 20 geothermal systems present at the surface, with a total estimated heat output of 4,200 MW (Bibby et al., 1995; Rowland & Simmons, 2012). 

## **2. Model Assumptions/Inputs** 

### **2.1. Heat Loss Timescale Calculation** 

We estimate the total heat loss (ΔHtotal) associated with cooling and crystallization of Mamaku magma as 

Map of the Central Taupō Volcanic Zone (TVZ) showing the where mMamaku is the mass of the Mamaku magma body (g), _h_  _h_ MELTS / _m_ MELTS, where h is the enthalpy of the system (J/g), hMELTS (J), and mMELTS (g) are the enthalpy values calculated in rhyolite-MELTS and the mass of the system used in rhyolite-MELTS, respectively (see Methods), and i is either 0% or qtz, which are the crystal fractions where the magma begins to crystallize and where quartz begins to crystallize, respectively. The quantity mMamaku is calculated by multiplying the density of Mamaku pumice (2.34 g/cm<sup>3</sup> ) by the estimated volume of the Mamaku ignimbrite (150 km<sup>3</sup> ). ΔHtotal is the change in enthalpy (in J) of the system from the starting temperature to the temperature reached when the desired crystal content of 10 wt% crystals, the maximum crystal content observed for Mamaku pumice, is reached. The temperature at the liquidus under fluid-saturated conditions calculated in the pure water case, which is when fluid is comprised of only H2O (hereby referred to as the wet liquidus), and the temperature when quartz first crystallizes, were both used as starting temperatures to obtain heat loss values for the entire history of the crystallizing melt and for the crystallization history of quartz. We assume that plagioclase is the phase present at the wet liquidus in our modelling. We consider the temperature at which quartz starts crystallizing because that leads to timescales directly comparable to Ti-in-quartz timescales. However, rhyolitic magmas are expected to saturate in at least one feldspar prior to quartz saturation (Gualda & Ghiorso, 2013a)—we use rhyolite-MELTS to estimate the total heat loss, and thus total crystallization time from the onset of feldspar crystallization. Figure 2a shows heat loss in rhyolite-MELTS for an average Mamaku ignimbrite composition (from Gualda et al., 2018), and Figure 2b shows heat loss for a range of Mamaku whole-pumice compositions, all of which are listed in Table 1. We note that a couple of compositions do not saturate in quartz prior to reaching 10 wt% crystals—these results can be explained by uncertainties in whole-rock and glass compositions, as well as consequent uncertainties in crystallization pressure. We contend that our calculations provide reasonable estimates for the total heat loss and total crystallization times for Mamaku magmas. 

**Figure 1.** Map of the Central Taupō Volcanic Zone (TVZ) showing the outline of the Rotorua caldera, which was active during the TVZ flare up from ∼350 to 240 ka. Outlines show the present-day extent of these flare up ignimbrites with the focus of this study, the Mamaku ignimbrite in green. Other flare-up units shown in gray. Information and inspiration for figure from Gravley et al. (2016). 

Once the total heat lost is calculated, the corresponding time is calculated using the following equation: 

where 

Q  q  Amagma (3) 

3 of 17 

KELLY ET AL. 

**Geochemistry, Geophysics, Geosystems** 

10.1029/2021GC009794 

<!-- Start of picture text -->
a<br>b<br><!-- End of picture text -->

**Figure 2.** Change in enthalpy per mass versus wt% crystals for (a) average Mamaku whole-pumice composition (at various pressures) and (b) individual pumice clasts from Mamaku whole-pumice compositions (Table 1). Changes in slope indicate the crystal content at which quartz begins to crystallize, referenced as qtz in the text. Enthalpy values are calculated during an isobaric cooling sequence performed using rhyoliteMELTS. All compositions are modeled at their own respective pressures at water oversaturation (i.e., enough water is added to maintain saturation throughout the simulation). Compositions in (b) without the change in slope did not crystallize quartz before 10 wt% crystals was reached, and only have a crystallization time calculated from, 0% to 10%. Compositions that yielded quartz crystallization have less enthalpy change over the entire crystallization history of the magma body. 

where q is the heat flux calculated using different models in the following sections. In the case of the Mamaku ignimbrite, we assume Amagma to be 150 km<sup>2</sup> , which corresponds to the surface expression of a circular magma body of ∼14 km diameter and a thickness of 1 km, consistent with the shape and area of the inferred Rotorua caldera (Milner et al., 2002). Crystallization in rhyolite-MELTS does not begin until the wet liquidus is reached, and quartz is typically not the first mineral phase to crystallize (see Gualda & Ghiorso, 2013a); therefore, two versions of tcool were calculated. The first considers the time to lose heat from the wet liquidus to 10 wt% crystals, and the other considers the time to lose heat from the onset of quartz crystallization to 10 wt% crystals (Gualda et al., 2018; Gualda, Pamukcu, et al., 2012; Pamukcu, Gualda, Ghiorso, et al., 2015). We make the distinction between wet liquidus and quartz crystallization times in order to compare calculated heat loss timescales from quartz saturation temperature to the aforementioned Ti-in-quartz chronometry timescales from Gualda et al. (2018). 

Our modeling approach is based on the research of Gualda et al. (2018, 2019), which suggest that Mamaku magmas were stored at relatively constant pressure, with crystallization due to heat loss over short timescales. Geobarometry results in Gualda et al. (2018, 2019) would be much different if crystallization primarily took place during ascent. Pumice textures and mineral zoning would also reflect this behavior. 

Section 3 of the paper aims to assess the viability of conductive and convective heat transfer as possible mechanisms for decadal to centennial heat transfer from a magma body by matching our calculated heat loss timescales, tcool, to the quartz crystallization timescales obtained by Gualda et al. (2018). Even if the timescales of Gualda et al. (2018) turn out to be inaccurate, our modeling constrains the mechanisms needed to cool magma on very short timescales. Our results also demonstrate the dramatic effects of convection in hydrothermal systems. 

Using the enthalpy change calculated with rhyolite-MELTS (ranging from ∼14 to 22 J/g) and the total estimated volume of the Mamaku magma (150 km<sup>3</sup> ) (Gravley et al., 2007), we find that the total heat loss needed to crystallize from wet liquidus to 10 wt% crystals is on the order of 3.4 × 10<sup>7</sup> –5.3 × 10<sup>7</sup> kJ/m<sup>2</sup> . In order for the system to lose the heat required to crystallize magma over 100 years, a heat flux of ∼11–17 W/m<sup>2</sup> is required, distributed over a surface area of 150 km<sup>2</sup> (Amagma). In order for the system to lose the required heat over 10 years, heat flux must be ∼110–170 W/m<sup>2</sup> . These calculated heat fluxes needed to crystallize the inferred Mamaku magma body from wet liquidus to 10 wt% crystals are subsequently referred to as qMELTS, with a range from 11 to 170 W/m<sup>2</sup> . 

As seen in Figure 2b, the MELTS modeling shows that not all modeled compositions crystallized quartz within 10 wt% crystals. This could mean that some samples do not have quartz, or, since crystal content is so low, it may be the case that there is quartz present in the sample but the modeling does not predict quartz crystallization. It could also be due to uncertainties in the whole-rock and glass compositions, as well as in the modeled pressure. 

For our study, we assumed that magma has been separated from the crystal mush and is being stored at shallow levels in the crust (Bégué et al., 2014; Gualda et al., 2018). In the TVZ, the mush is presumed to sit at and below the brittle-ductile transition (BDT) (Barker et al., 2021; Rowland & Simmons, 2012). The magmas that erupted to form the Mamaku ignimbrite are stored above the mush after extraction, demonstrated in Gualda et al. (2019). Therefore, we assume cooling of the extracted magma from the crystal mush and 

4 of 17 

KELLY ET AL. 

**Geochemistry, Geophysics, Geosystems** 

10.1029/2021GC009794 

#### **Table 1** 

|_Whole-Pumice Composit_|_ions of Mama_|_ku Ignimbr_|_ite Samples_||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|
||Average<br>Mamaku<sup>a</sup>|GL2306-<br>B<sup>b</sup>|GL2306-<br>C<sup>b</sup>|GL2306-<br>D<sup>b</sup>|GL2306-<br>E<sup>b</sup>|GL2306-<br>F<sup>b</sup>|GL2306-<br>G<sup>b</sup>|GL2306-<br>H<sup>b</sup>|GL2306-<br>K<sup>b</sup>|GL2306-<br>M<sup>b</sup>|GL2306-<br>N<sup>b</sup>|GL2306-<br>O<sup>b</sup>|
|SiO2|76.90|74.65|73.32|73.30|74.38|70.19|74.62|73.83|75.05|72.88|71.50|74.84|
|TiO2|0.15|0.18|0.20|0.20|0.18|0.37|0.18|0.20|0.17|0.22|0.27|0.19|
|Al2O3|13.10|15.63|16.93|17.11|15.89|17.76|15.55|16.15|15.18|17.40|17.63|15.15|
|Fe2O3|1.41|0.16|0.28|0.29|0.26|0.45|0.25|0.28|0.24|0.30|0.36|0.28|
|FeO<sup>g</sup>|0.01|1.31|1.33|1.38|1.19|2.12|1.16|1.28|1.11|1.45|1.68|1.26|
|MnO|0.04||||||||||||
|MgO|0.10|0.16|0.19|0.17|0.16|0.41|0.15|0.18|0.16|0.20|0.27|0.18|
|CaO|0.73|0.81|0.92|0.89|0.81|1.50|0.86|0.91|0.82|0.91|1.13|0.92|
|Na2O|3.96|3.76|3.66|3.70|3.79|4.10|3.88|3.88|3.82|3.74|3.78|4.44|
|K2O|3.73|3.13|2.94|2.75|3.15|2.75|3.15|3.06|3.25|2.69|3.12|2.54|
|P2O5|0.02||||||||||||
|H2O|4.89|4.85|5.21|5.09|5.17|4.49|4.97|4.56|5.20|5.14|6.17|4.89|
|_P_(MPa)<sup>c</sup>|50–150|144|160|152|160|124|149|128|163|154|211|138|
|_l_(km)<sup>d</sup>|∼1.8–5.5|∼5.3|∼5.9|∼5.6|∼5.9|∼4.6|∼5.5|∼4.7|∼6|∼5.7|∼7.8|∼5.1|
|Δ_H_from quartz<br>saturation (J/g)<sup>e</sup>|∼15–23|7.05|||7.13||3.64||7.51||||
|Δ_H_from wet liquidus<br>(J/g)<sup>f</sup>|∼30–49|62.75|75.57|73.63|62.35|72.90|63.35|72.37|52.28|75.45|69.45|66.24|

aAverage composition is taken from Gualda et al. (2018). bAll samples labeled GL2306 are individual pumice clasts analyzed as part of this study. cStorage pressure of magma body used in rhyolite-MELTS modeling.<sup>d</sup> Corresponding depth of magma body calculated from _P_ using _P_ = _ρgl_ .<sup>e</sup> Calculated change in enthalpy from rhyolite-MELTS from quartz saturation until 10 wt% crystals.<sup>f</sup> Calculated change in enthalpy from wet liquidus temperature until 10 wt% crystals. gInitial values for FeO and Fe2O3 were determined using rhyolite-MELTS from measured total FeO values by assuming _fO2_ equal to that of the Ni-NiO buffer. 

not magma embedded within, which would require cooling of the mush, increasing cooling timescales as a result. 

### **2.2. Current Surface Heat Flux in the TVZ** 

Previous estimates of the overall heat budget of the TVZ (Hochstein, 1995) show that the present-day heat output is between 3,000 and 4,000 MW across the entire rift, which is approximately 200 km long and 25 km wide. With these values, the rift-average heat flux is between 0.6 and 0.8 W/m<sup>2</sup> . Alternative estimates (Dempsey, 2012; Weir, 2009) take into account the average area of a geothermal field, the average power of a geothermal field, and how many geothermal fields fit into the footprint of a magma body. The average heat flux through a geothermal field is calculated as 1.4 W/m<sup>2</sup> using these assumptions (Dempsey, 2012; Weir, 2009). Using 0.7 and 1.4 W/m<sup>2</sup> as the lower and upper bounds, respectively, for our calculations, subsequently referred to as qTVZ, calculated timescales of heat loss, tcool, are on the order of millennia, which are too slow compared to the crystallization timescales from quartz geochronometry from Gualda et al. (2018) (see above). 

## **3. Testing Viability of Models With Storage Timescales as Determined by Quartz Crystallization** 

### **3.1. Pure Conduction** 

For this model, we assume heat conduction occurs through a homogeneous layer of crustal material between the Mamaku magma body and the surface of the TVZ (Figure 3a). Heat is lost over a static, linear geothermal gradient according to Fourier's Law. In this case, the flux of heat per unit area is constant and given by Carslaw and Jaeger (1984) 

5 of 17 

KELLY ET AL. 

**Geochemistry, Geophysics, Geosystems** 

10.1029/2021GC009794 

**Figure 3.** Conceptual models for analytical heat flow solutions used in this study. (a) Schematic diagram of constant conduction based on Equation 4 from Carslaw and Jaeger (1984). Depths are in line with those reflected in Figure 2a. (b) Schematic diagram of combination of conductive and convective heat flow, with conductive layer directly above the magma body overlain by a hydrothermal component first detailed in Section 3.2 and expressed by Equation 7. The top of the conductive layer is constrained to 300–350°C (Dempsey et al., 2018; Rowland & Simmons, 2012). Depths reflected are constrained by distribution created from GL2306-G glass average and standard deviation, and are the depths used in Monte Carlo modeling of Equations 7 and 8. To see more details, see Section 6.3. (c) Schematic diagram of Rainfall-influenced hydrothermal model (Equation 5) (Dempsey et al., 2018; Weir, 2009). All models output a heat flux in W/m<sup>2</sup> , and timescales are calculated as outlined in the text. 

where  is the thermal conductivity (W m<sup>−1</sup> K<sup>−1</sup> ) of the crust, T1, T2 are the temperatures of the top surface and the magma body (K), and  is the depth to the top of the magma body (in m). An average value for  of 2.3 W m<sup>−1</sup> K<sup>−1</sup> is taken from Eppelbaum et al. (2014). The values for T2 and T1 are respectively the wet liquidus, which ranges between 781 and 861°C, and 25°C. The temperature of the magma body is taken to be constant (in practice, it varies about 2%—up to 15 K—over the range 0–10 wt% crystals; see Gualda et al., 2018). Thickness  corresponds to the depth of the magma body below the surface, which is obtained from crystallization pressures from the rhyolite-MELTS geobarometer (Gualda & Ghiorso, 2014; Gualda et al., 2018) using the equation P  gl, where P is the storage pressure of the magma body prior to eruption (see Tables 1 and 2),  is density of 2.75 × 10<sup>3</sup> kg/m<sup>3</sup> , and g is 9.8 m/s<sup>2</sup> (e.g., 100 MPa corresponds to ∼3.7 km depth). 

Using this model, timescales obtained are millennial in magnitude, with all timescales from quartz saturation temperature over 1,000 years (Figure 4a), consistent with the calculations of Gualda, Pamukcu, et al. (2012). The effect of variable storage pressure is negligible over the relatively narrow range of pressures estimated for the Mamaku magma body in the shallow crust (Bégué et al., 2014; Gualda et al., 2018), with timescales over 1,000 years for all pressures modeled. 

These timescales are not consistent with the primarily decadal timescales of quartz crystallization from Gualda et al. (2018), indicating that conduction is not the dominant heat transfer mechanism during this period. The heat flux, q, for the typical storage depths shown in Table 2 and calculated by conduction ranges 

6 of 17 

KELLY ET AL. 

**Geochemistry, Geophysics, Geosystems** 

10.1029/2021GC009794 

||_σ_for<br>GL2306-G<br>Glassb|0.14|0.04|0.14|0.15|0.02|0.04|0.05|0.06||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|
||GL2306-O|77.66|0.12|12.65|0.59|0.06|0.46|3.90|4.57|10|138|∼5.1|
||GL2306-N|76.94|0.18|13.01|0.87|0.07|0.67|4.15|4.11|10|211|∼7.8|
||GL2306-M|77.57|0.13|12.73|0.51|0.02|0.39|4.04|4.60|10|154|∼5.7|
||GL2306-K|77.41|0.11|12.74|0.71|0.03|0.43|4.09|4.47|10|163|∼6|
||GL2306-H|77.64|0.15|12.61|0.71|0.03|0.44|3.95|4.47|10|128|∼4.7|
||GL2306-Ga|77.55|0.14|12.75|0.62|0.02|0.47|3.99|4.46|10|149|∼5.5|
||GL2306-F|77.60|0.14|12.59|0.83|0.05|0.54|3.87|4.37|10|124|∼4.6|
|_umice Samples_|GL2306-E|77.50|0.13|12.75|0.66|0.03|0.44|4.09|4.41|10|160|∼5.9|
|_aku Ignimbrite P_|GL2306-D|77.53|0.11|12.79|0.56|0.02|0.39|4.04|4.55|10|152|∼5.6|
|_ositions of Mam_|GL2306-C|77.55|0.12|12.75|0.47|0.03|0.47|3.90|4.71|10|160|∼5.9|
|_trix Glass Comp_|GL2306-B|77.50|0.12|12.67|0.78|0.04|0.42|4.04|4.44|10|144|∼5.3|
|**Table 2**<br>_Average Ma_||SiO2|TiO2|Al2O3|FeO|MgO|CaO|Na2O|K2O|H2O|_P_(MPa)c|_l_(km)d|

from 0.3 to 1.1 W/m<sup>2</sup> , significantly smaller than qMELTS, and similar in magnitude to qTVZ. The calculation can also be inverted: if observations of qTVZ are used to compute equivalent storage depths for pure conduction heat loss, these range between 1.3 and 3 km. This is much shallower than storage depths given in Table 2. Storage depths of magma bodies that crystallize by pure conduction in decades or centuries (qMELTS) range between 13 and 130 m, which are geologically unreasonable for large, silicic magma bodies. 

### **3.2. Convection** 

The inability of a pure conductive layer to lose enough heat on the timescales inferred from geochronometry indicates that the appropriate heat extraction mechanism likely incorporates a hydrothermal convective component. This is particularly likely in the case of the TVZ (but also many other volcanically active areas), where geothermal systems are ubiquitous and well developed (Bibby et al., 1995; Rowland & Simmons, 2012). In such a convective system, conduction occurs in a lower boundary layer, the viscoelastic shell, above the magma body until the boiling temperature of the hydrothermal fluids is reached at ∼300–350°C (Dempsey et al., 2018; Weir, 2009) (Figure 3b). Above this layer, advection of hot fluids leads to a region in which temperature is nearly constant and convective heat transfer is dominant. Close to the Earth's surface, temperature rapidly drops to a surface atmospheric value (Rowland & Simmons, 2012). Introducing fluid advection allows for heat extraction systems that are more efficient than conduction alone. 

Energy conservation requires that, in steady state, each of the layers described above (lower conductive layer, middle advective layer, and upper boundary layer) transfer heat at the same rate. For the lower conductive layer, temperature boundary conditions are well constrained, at ∼750°C in contact with the magma body (Gualda et al., 2018) and at 300–350°C at the BDT temperature of quartzo-feldspathic crust, after the onset of quartz plasticity at 300°C (Scholz, 1988), assumed to be the lower permeable limit accessible by hydrothermal fluids (Dempsey et al., 2018; Rowland & Simmons, 2012). There is also the work of Fournier (1999), in which similar temperatures of the BDT around 360 °C are estimated, where hydrostatic pressure is dominant and strain rates are in line with common values of actively deforming areas. However, Fournier (1999) also shows that times where strain rate increases contribute to a hotter BDT of 870°C, but these times of higher strain would occur during dome emplacement or collapse due to eruption, processes not associated with storage. In order to crystallize 10 wt% crystals in the Mamaku magma on decades to century timescales, the lower conductive layer would need to be ∼5.9–102 m thick. This is (mostly) an order of magnitude smaller than the estimate of Dempsey et al. (2018), who calculated a value of ∼450 m. Depending on the temperature at the brittle ductile transition, a minimum magma volume of ∼100 km<sup>3</sup> can be cooled over 100 years with a lower conductive layer on the order of 100 m thick. A conductive layer thickness approaching ∼450 m would be sufficient for magma of ∼22–25 km<sup>3</sup> . This implies that the magma system may likely be divided between different lenses that are smaller in volume than the modeled 150 km<sup>3</sup> , with smaller volumes cooling over shorter timescales than larger volumes that are cooled at the same rate. 

7 of 17 

KELLY ET AL. 

**Geochemistry, Geophysics, Geosystems** 

10.1029/2021GC009794 

<!-- Start of picture text -->
a<br>b<br><!-- End of picture text -->

<!-- Start of picture text -->
b<br><!-- End of picture text -->

**Figure 4.** Timescales of heat loss for the pure conduction and rainfall convection analytical models. (a) Results of conductive heat flow modeling using Equation 4 from Carslaw and Jaeger (1984). (b) Results of Rainfall modeling using Equation 5 from Weir (2009) and Dempsey et al. (2018). All timescales were calculated from the enthalpy changes of the Average Mamaku compositions in Figure 2a. The wet liquidus line shows the time to lose heat over the entire crystallization history from 0 to 10 wt% crystals of the magma body, whereas the quartz line shows the time to lose heat from quartz saturation temperature until 10 wt% crystals. 

### **3.3. Rainfall Convection** 

It is estimated that >85% of modern convecting hydrothermal fluids in the TVZ are meteoric rather than magmatic in origin (Giggenbach, 1995). As such, any analysis of heat loss in the TVZ would be remiss not to include meteoric fluid-driven hydrothermal heat flow. Weir (2009) puts forth a “once-through” model of the infiltration of cold rainwater and upwelling of heated water into a TVZ geothermal system, where the rainwater only passes through the system once. Circulation is controlled by density differences between the downwelling meteoric water and upwelling geothermal fluids, as well as how much water is able to infiltrate the system. This model was expanded upon by Dempsey et al. (2018), who introduced the idea that the downwelling region of a geothermal system is 7 times larger than the upwelling region (Dempsey et al., 2012), when considering that upwelling regions in the TVZ range between 5 and 35 km<sup>2</sup> (Bibby et al., 1995). They derived an expression for the heat flux 

where hu is the upflow enthalpy of ∼1,550 kJ kg<sup>−1</sup> for liquid water at 350°C at 4–6 km depth, d is the density of the downwelling fluid (1,000 kg m<sup>−3</sup> ), R is the regional rainfall rate (∼1,200 mm yr−1), 7 / 8 is the ratio of geothermal fields an area can support, and  is the fraction of rainfall that infiltrates the groundwater in the TVZ. Based on geothermal circulation modeling, Kissling and Weir (2005) and O'Sullivan et al. (2009) have determined  to be in the range 5%–10% (see Figure 3c). For simplicity, the yearly rainfall rate and fraction of rainfall infiltration are assumed to be constant in our treatment. 

Timescales of heat loss from quartz saturation temperature obtained using this model are between ∼200 and 300 years, and timescales of heat loss from wet liquidus temperature are between 400 and 700 years (Figure 4b). These timescales are similar in order of magnitude to some timescales obtained by Ti-in-quartz geochronometry (Gualda et al., 2018), but, again, are longer than the majority of quartz crystallization timescales. Heat loss timescales over the entire lifetime of the magma body, however, reach centennial timescales, which is within the estimated range for the estimated total residence time of a crystal-poor, rhyolite magma body in the TVZ (Gualda et al., 2018). Variations in timescales due to variable storage pressure are, once again, negligible for the range of shallow pressures used. 

There is uncertainty in the R<sup></sup> today as well as in the past. Results above are calculated using   0.1, or 10%. Increased values of  on the order of 20%–25% would lead to heat loss on decadal timescales; however, it is not clear whether such values could be reasonably attained. A present-day estimate for  is between 5% and 10% (Kissling & Weir, 2005; O'Sullivan et al., 2009), and any changes are difficult to quantify and would involve taking into account the tectonic state of the TVZ, which is beyond the scope of this paper. This approach also assumes that the state of permeability in the crust is the same as it was prior to eruption of the Mamaku. It is possible that a period of accelerated rifting in the upper crust (cf. Gravley et al., 2016) could enhance extensional fault/fracture controlled permeability (see Rowland & Sibson, 2004) relative to the background tectonic state, leading to an increase in . This rainfall-controlled model also brings up the question of the importance of the climate in the region where the magma body is located. If there is a high annual rainfall rate, heat flux would be higher and tcool would be shorter. If this model is chosen, then one must consider the impact of the climate in the region of interest in any results. 

8 of 17 

KELLY ET AL. 

**Geochemistry, Geophysics, Geosystems** 

10.1029/2021GC009794 

**Table 3** 

_List and Description of Parameters Used in Monte Carlo Modeling_ 

|Parameter|Range|Units|Distribution|Description|
|---|---|---|---|---|
|_Vz_|8.7 × 10<sup>−10</sup>–4.8 × 10<sup>−9</sup>|m/s|Uniform|Darcy velocity of groundwater|
|_β_|−10 to −50|dimensionless||Peclet number|
|_κ_|2.0–2.5|W/mK|Uniform|Thermal conductivity|
|_c_0|4.2 × 10<sup>3</sup>|J/K||Specific heat of groundwater|
|_ρ_0|800–900|kg/m<sup>3</sup>|Uniform|Density of groundwater|
|_L_|∼4–7|km|Normal|Storage depth of magma body|

_Note_ . All values from Dempsey (2012), except _κ_ , which is from Eppelbaum et al. (2014). Unknown distributions are assumed uniform according to the Principle of Indifference (Keynes, 1921). The distribution of _L_ is normal to due the Gaussian nature of the distribution of the GL2306 glass compositions, whose average and standard deviation are shown in Table 2. 

### **3.4. Groundwater Flow** 

One key shortcoming of the convective model outlined above is that it does not take into account how the speed of the groundwater moving through the crust impacts heat flow seen at the surface. This model characterizes the hydrothermal convective component of the system illustrated in Figure 4b. The model characterizes heat flux as a derivation of the one-dimensional advection-diffusion equation (Ardid et al., 2021; Bredehoeft & Papadapulos, 1965), where 

Heat loss is assumed to occur in a steady state (not time dependent) in the vertical ( _z_ ) direction. Heat is assumed to be constant in the horizontal directions. The final characterization of the flux is 

and is controlled by a Peclet number, , where 

which is a description of the ratio of advective to diffusive transport that is influenced by the Darcy velocity of groundwater in the vertical direction, v . The temperatures z TL and T0 represent the temperature of the convecting fluids at depth (300–350°C) and temperature of the surface (25°C), respectively. For a full derivation, see Methods in Supporting Information S1. Values for  in the TVZ range from −10 to −50 (Dempsey, 2012). All values used in the calculation of  are presented in Table 3. Heat flux is assumed to be averaged over the total surface area, so the solution is one-dimensional and only concerned with the change in flux in the vertical direction. Due to the uncertainty in calculating the storage pressure (parameter L from above) of the magma body a Monte Carlo analysis of storage pressures was performed (see Methods, Section 6.3), where a distribution of storage pressures of magma compositions is calculated using the rhyolite-MELTS geobarometer. There are also uncertainties in the other parameters used to calculate the Peclet number , so a Monte Carlo analysis of the parameter was also performed to generate a distribution of heat loss timescales calculated with the groundwater flow model (see Methods, Section 6.3). A random number was drawn from the distributions of each parameter 10,000 times and used to calculate , which was then used to calculate a corresponding heat flux and heat loss timescale tcool. 

From the groundwater flow model we calculate heat loss timescales from quartz saturation to 10 wt% crystals of 0–1,000 years. The minimum is 0 because for GL2306-M, quartz saturated at ∼10 wt% crystals for a water fraction relative to saturation of 0.1, so there was no enthalpy change. The next lowest timescale value is ∼18 years. The distribution of heat loss timescales calculated using the groundwater flow model from quartz saturation temperature until 10 wt% crystals shows values skewed toward decades to one century 

9 of 17 

KELLY ET AL. 

**Geochemistry, Geophysics, Geosystems** 

10.1029/2021GC009794 

<!-- Start of picture text -->
a<br>b<br><!-- End of picture text -->

**Figure 5.** Distribution of heat loss timescales from (a) quartz saturation and (b) from wet liquidus temperature generated by the Monte Carlo method of the groundwater heat flow model. 

(Figure 5a). The distributions of timescales from wet liquidus temperature to 10 wt% crystals are also skewed, being skewed toward fewer than multiple centuries (Figure 5b). 

The skewed heat loss timescales are not captured by the rainfall and conduction solutions put forward in the previous sections, and the asymmetrical behavior of the distribution makes qualitative sense because time cannot be negative, so a concentration of values around zero is not surprising. In this sense, these results show the behavior of the system when taking into account that the distributions of the parameters used to calculate the timescales are Gaussian and uniform. The results indicate that the groundwater flow model can be utilized for further understanding heat loss of a magma body on decadal timescales. 

## **4. Discussion** 

### **4.1. Conduction Versus Convection** 

The initial goal of this study was to pursue conductive and convective heat flow models to determine which mechanism is more likely to contribute to heat loss from a magma body on timescales on the order of decades to centuries. When comparing conductive and convective heat flow, the results indicate that decadal to centennial-scale heat loss cannot happen as a result of pure conductive heat flow over the entire crustal thickness above a large magma body, with all timescales for heat loss by conduction leading to crystallization times greater than 1,000 years. The heat flux of the rainfall convection model is much closer to qMELTS than heat flux from the constant conductive model (∼5 W/m<sup>2</sup> vs. 1.3 W/m<sup>2</sup> ). Therefore, a variant of a convective system characterized by conductive heat transfer in the crust around the magma body but dominated by convective heat transfer through hydrothermal fluids is more likely to be the mechanism by which heat was extracted from the magma body that led to the caldera-forming eruption of ∼150 km<sup>3</sup> of Mamaku magma, in order to be consistent with crystallization timescales derived from quartz crystals. We further explore—in addition to the rainfall convective model—a groundwater flow model, which is consistent with the convective scenar- 

io proposed above, but which allows us to better investigate some of the variations possible in convective heat loss. In an area like the TVZ, with ∼20 hydrothermal systems (Bibby et al., 1995; Rowland & Simmons, 2012), this heat extraction mechanism is likely to have an effect on relatively shallow silicic magma bodies cooling in the crust. 

Variations in water content and magma composition can lead to variations in the value of ΔHtotal. A particularly striking aspect of the behavior of convective heat loss timescales using the groundwater flow model is that there is persistence in the order of magnitude of decadal to centennial values for tcool (Figure 6). Even when crystallization in rhyolite-MELTS is extended until 20 wt% crystals, which leads to an increase in ΔHtotal, the majority of values of tcool calculated using flux from Equation 7 are centennial in magnitude (Figure 6c). 

In order to analyze the probability that calculated heat loss timescales from the groundwater-controlled hydrothermal model exceed 100 years, we show the distribution of the exceedance probabilities of heat loss times, tcool, whose enthalpy changes are calculated from quartz saturation qtz until 10 wt% crystals. Figure 7a demonstrates that 50% of timescales are less than 100 years, and 90% are less than 280 years. These timescales are in line with the evidence of short crystallization and storage on the order of decades to centuries from Gualda et al. (2018), drawing an important connection between a crystallizing magma body and the crustal processes that remove heat from it. 

10 of 17 

KELLY ET AL. 

**Geochemistry, Geophysics, Geosystems** 

10.1029/2021GC009794 

<!-- Start of picture text -->
a<br>b<br>c<br><!-- End of picture text -->

**Figure 6.** Average timescales of heat loss from quartz saturation temperature. Timescales were obtained using the groundwater heat flow model (Equation 7) for a range of compositions and water fractions. Timescales represent crystallization from quartz saturation until (a) 10, (b) 15, and (c) 20 wt% crystals. Different crystal contents were considered to allow for crystallization of quartz in all compositions. Colors of contours correspond to the timescale of heat loss from quartz saturation averaged over all calculations from the Monte Carlo analysis. 

### **4.2. Coupling of Magmatic and Hydrothermal Heat Loss** 

A question that arises from this research is how the heat is distributed throughout the magma body. If it is assumed that the entire magma body is well-mixed, as we do, then, based on the maximum value of qTVZ of 1.4 W/m<sup>2</sup> (Dempsey, 2012; Weir, 2009), the thickness of a magma body that loses the required amount of heat in 100 years is ∼100 m, and the thickness of a magma body that loses the required heat in 10 years is ∼10 m. Both values give a magma thickness 1–2 orders of magnitude less than the assumed 1 km thickness, seemingly indicating that a wellmixed magma body of volume 150 km<sup>3</sup> could not sufficiently cool in the required time. However, the actual subsurface expression of the Mamaku magma body could be a patchwork of multiple smaller magma bodies (Bégué et al., 2014), with different magma lenses or different volumes that are not arranged in any particular way but are within the footprint of the caldera structure, rather than the simpler assumption of a single pancake-shaped magma body (G. F. Cooper et al., 2012; see also Cashman & Giordano, 2014; Gualda & Ghiorso, 2013b). This patchwork arrangement, where smaller magma bodies of sizes such as 30 km<sup>3</sup> could possibly allow for these smaller magma bodies to cool, erupt, and possibly trigger eruption of other magma bodies in the system. This could resolve the disagreement between our calculations of the size of the lower conductive viscoelastic shell above the magma body and the calculations of Dempsey et al. (2018). These magma bodies could be spread out over an area greater than 150 km<sup>3</sup> , which would increase total heat output if constant heat flux is assumed; or they could span a thickness greater than 1 km, which could impact the depth extent incorporated into our calculations. If the magma is believed to be one contiguous body, the effect of transient permeability on decreasing the size of the viscoelastic shell could be investigated (Mittal & Richards, 2019). 

### **4.3. Is There a Mamaku-Sized Caldera-Forming Magma Body Beneath Rotorua Now?** 

Comparison of heat flow timescales to quartz crystallization timescales shows that hydrothermal systems are the only mechanism able to remove heat from a crystallizing magma body at sufficiently fast rates. We now explore more closely the effect of Monte Carlo analysis of the groundwater heat flow model on the calculated heat output, Q, using values calculated from the current TVZ as constraints for our analysis. Figure 7b shows the distribution of heat output calculated as a result of the Monte Carlo analysis of the groundwater hydrothermal model that yield heat loss timescales of 100 years or less. The mean and median heat output values are ∼1,000 MW. When looking at the distribution in Figure 7b, the values for heat output range from ∼400 to 1,600 MW, or 1 to 4 times the heat output from the current Rotorua hydrothermal system. However, the majority of power output values are 2–3 times the size of the current Rotorua system (470 MW) (Bibby et al., 1995), which indicates it is not very likely that a Mamaku-sized magma body is currently present in the crust and cooling on the order of decades below the Rotorua caldera. It is likely that larger heat output values could be typical at times when a large, silicic magma body is being stored in the crust. 

The Rotorua caldera is currently defined by two ellipses that create a structural boundary that is 21 km by 20 km, or 420 km<sup>2</sup> (Milner et al., 2002). If 

11 of 17 

KELLY ET AL. 

**Geochemistry, Geophysics, Geosystems** 

10.1029/2021GC009794 

<!-- Start of picture text -->
a<br>b<br><!-- End of picture text -->

**Figure 7.** Distribution of timescales and heat output from Equation 7. (a) Cumulative distribution (blue) and exceedance probability (orange) functions of timescales of heat loss from quartz saturation temperatures. Timescales used were calculated using the groundwater heat flow model (Equation 7) and shown in histogram form in Figure 5a (see above). All times that were 0 were not incorporated into the function calculations. Virtually all timescales are less than 1,000 years, with 90% less than 280 years. (b) Distribution of heat outputs that contribute to heat loss over timescales less than 100 years from quartz saturation temperature until 10 wt% crystals (from Figure 5a). Mean and median values of heat output are around 1,000 MW. Heat outputs are on the order of 1–4 times the size of the current Rotorua hydrothermal system. Heat outputs shown were calculated using the groundwater heat flow model (Equation 7), by assuming flux over a surface area of 150 km<sup>2</sup> . 

we assume that this 470 MW is evenly spread over the entire area of the caldera, according to Bibby et al. (1995), then the heat flux over the area of the Rotorua caldera would be 1.12 W/m<sup>2</sup> , more in line with calculated conductive heat fluxes and qTVZ. If we assume that heat is being lost over the entire extent of the Rotorua caldera (420 km<sup>2</sup> ), at a rate of 1.12 W/m<sup>2</sup> , then calculated times of heat loss are around 500 years from quartz saturation, and up to 1,000 years for the entire crystallization history of the magma body. This indicates that, in order to cool a magma body on the order of decades to centuries using the current Rotorua caldera area and total heat output of hydrothermal systems within the caldera, a current magma body should be 30 km<sup>3</sup> or less in size, 1/5 of the assumed 150 km<sup>3</sup> magma body volume. 

### **4.4. Implications for the Current TVZ and Beyond** 

While heat flow due to hydrothermal convection appears to be essential for decadal heat loss timescales from a magma body the size of that inferred for eruption and deposition of the Mamaku ignimbrite, it is important to note that the presence of hydrothermal systems of that magnitude may not mean that there is a large, eruptible, silicic magma body being stored in the shallow crust. Currently, there are areas within the TVZ that exhibit heat output similar to the mean heat output from the groundwater heat flow model displayed in Figure 7b. The Waimangu, Waiotapu, and Reporoa hydrothermal systems, which are located roughly 20 km southeast of the current Rotorua system, together have a heat output of ∼1,000 MW over an area of ∼110 km<sup>2</sup> (Bibby et al., 1994, 1995; Rowland & Simmons, 2012). Rowland and Simmons (2012) indicate that heat flow in the TVZ is relatively stable despite frequent volcanic eruptions. They also champion the idea of a “hot zone” in the lower to middle crust, where basaltic magmas intrude and contribute to the generation and eruption of silicic mush that drives the present-day heat flow in the region (Rowland & Simmons, 2012). Therefore, heat flow values similar to those inferred to have existed above the Mamaku magma body just prior to eruption do not necessarily indicate shallow magma storage. 

An understanding of when a large, silicic magma body is being stored in the crust is essential for distinguishing between relatively benign hydrothermal activity and activity that indicates the possibility of a large, caldera-forming eruption. Whether that be through seismicity, ground deformation, anomalous gas signatures in hydrothermal systems that indicate new magma, possible elevated heat flow measurements, or the generation of new hydrothermal systems, accumulation of a magma body of ∼150 km<sup>3</sup> will likely bring about change that can be viewed at the surface. If we can characterize and quantify these changes in the hydrothermal system, we may have the potential to detect a large volume of silicic magma being stored in the shallow crust and gauge its eruptibility with respect to time. 

## **5. Conclusion** 

We have compared the heat loss timescales of a magma body under the influence of conductive and convective heat flow models. Convective models that are governed by hydrothermal systems give heat loss timescales closest to crystallization times from quartz geochronometry, with a heat output on the order of 1,000 MW required for removing heat from a large, silicic magma body on the order of decades to centuries. This ephemeral nature of the Mamaku magma body holds regardless of storage pressure, water saturation, 

12 of 17 

KELLY ET AL. 

**Geochemistry, Geophysics, Geosystems** 

10.1029/2021GC009794 

or magma composition in our study, and is reflected in the fact that 50% of calculated timescales from the groundwater flow model are less than 100 years, and 90% are less than 280 years. 

The question of how heat is lost throughout the magma body brings up the possibility of multiple magma bodies of smaller size crystallizing on the order of decades, rather than a single well-mixed magma body 150 km<sup>3</sup> in volume. This can resolve discrepancies between thicknesses of the viscoelastic shell around the magma body, as well as decrease the heat output required of a hydrothermal system to lose heat on the order of decades. Future studies could explore this relationship further by coupling magma bodies and hydrothermal systems in more detail. 

The current heat output of the Rotorua hydrothermal system is 2–3 times smaller than the majority of power outputs calculated using the groundwater heat flow model that used the current variability in the modeling parameters in the TVZ, so the presence of a Mamaku-sized magma body cooling on the order of decades below Rotorua at present is not likely. There are three closely spaced hydrothermal bodies in the TVZ that collectively have a heat output on the order of 1,000 MW, but the driver of these systems is thought to be from a deeper zone of the crust than the depth of a shallowly stored silicic magma in the TVZ. Thus, the challenge is knowing when such large silicic magma bodies are being stored in the shallow crust so that we can apply our increasing understanding of their eruptible timescales. Through knowledge of distribution of magma bodies as single entities or a patchwork of bodies in the crust as well as when silicic magma is collecting in the crust, future studies could utilize the results of this investigation in relation to cooling of magma bodies of different size under different heat flow regimes in order to improve the ways we monitor restless calderas in New Zealand and elsewhere. 

## **6. Materials and Methods** 

### **6.1. Samples** 

Whole-rock compositions for 11 samples of Mamaku ignimbrite pumice were used as part of this study (Table 1). The samples were collected from an outcrop along Sutcliffe Rd. (−38.15788, 176.01186, WGS 84). 

Whole-pumice compositions were obtained using a Siemens SRS3000 dispersive X-ray Fluorescence Spectrometer with a 3 kW Rhodium X-ray tube at Verum Group Ltd (Spectrachem) in New Zealand, following their standard procedures. Samples were ground to a fine powder with a tungsten-carbide ring mill, dried at 110°C overnight at a minimum, ignited for 1 h at 1000°C, then combined with a proprietary mixture of 12:22 lithium tetraborate and lithium metaborate X-ray flux. The sample was mixed and heated to 1100°C, producing a glass melt used to form a flat disk. The glass was analyzed against an analytical calibration developed for alumino-silicate materials containing Fe2O3, MnO, TiO2, CaO, K2O, P2O5, SiO2, Al2O3, MgO, and Na2O, where each elemental analytical line is calibrated separately. Whole-pumice compositions are shown in Table 1. 

Glass compositions were obtained with a Tescan VEGA3 scanning electron microscope with an Oxford X-max 50-mm<sup>2</sup> energy-dispersive spectrometry (EDS) system for major-element compositional analysis. We used a 15 kV beam, 15 mm working distance, with specimen currents on the order of 3 nA for EDS analysis. We used live acquisition times of 25 s, so as to minimize Na loss but yield measurable amounts (i.e., above quantification limit) of minor elements, particularly Ti and Mg, which appear in concentrations between 0.1 and 0.5 oxide wt%. Acquisition areas correspond to rectangular or irregular regions typically ∼10 µm in linear dimension, drawn to avoid crystals or epoxy. We set process time to 4 on the AZtec acquisition and quantification software from Oxford Instruments, and O was calculated by stoichiometry, with quantification performed using internal standards as customary in EDS analysis. Each analysis presented is the average of 10–20 individual analyses from glass in each pumice clast. The quality of the resulting analysis was checked by comparing analyses performed under the same conditions listed above for a rhyolitic glass synthesized in the laboratory at Vanderbilt University using the U.S. Geological Survey reference powder RGM-1 (Gualda et al., 2018; see also Pamukcu, Gualda, Ghiorso, et al., 2015). Matrix glass compositions are shown in Table 2. 

Following Gualda et al. (2018), we use rhyolite-MELTS geobarometry (Gualda & Ghiorso, 2014; see also Bégué et al., 2014; Gualda et al., 2018), to calculate storage depths of magma bodies prior to eruption, which 

13 of 17 

KELLY ET AL. 

**Geochemistry, Geophysics, Geosystems** 

10.1029/2021GC009794 

is an important parameter for some of our models. Errors associated with pressure estimates are discussed by Gualda, Ghiorso, et al. (2012) and Pitcher et al. (2020). Resulting storage pressures are also shown in Table 2. 

### **6.2. Total Enthalpy Loss Modeling** 

In order to model the total heat loss of the Mamaku magmatic system over its cooling and crystallization history, we used rhyolite-MELTS, a set of thermodynamic models developed to model fluid-bearing silicic magma systems (Gualda & Ghiorso, 2014; Gualda, Ghiorso, et al., 2012; see also Asimow & Ghiorso, 1998; Ghiorso & Sack, 1995). Enthalpy loss as a function of crystallization is assumed to be a proxy for crystallization time of the magma body. Heat loss is assumed to be in a steady state, causing enthalpy to be correlated with time by a constant heat loss rate that is spatially averaged over the surface expression of the Mamaku magma body. Rhyolite-MELTS is used to calculate isobaric cooling sequences, which yield a complete description of the crystallization path (phases present, their compositions, thermodynamic properties of phases and system, etc.) for a given system composition for a set temperature interval, at a specified pressure. Whole-rock compositions of pumice are used because pumice is believed to be an erupted parcel of magma, recording the composition of the system prior to crystallization. 

To explore the effect of pressure on the total enthalpy loss, separate isobaric cooling sequences are calculated using the average Mamaku pumice composition of Gualda et al. (2018), at pressures from 50 to 150 MPa in 25 MPa pressure steps. The magma is assumed to be initially fluid saturated. Initial values for FeO and Fe2O3 were determined using rhyolite-MELTS from measured total FeO values by assuming _fO2_ equal to that of the Ni-NiO buffer. Calculations are allowed to continue until a maximum of 10 wt% crystals is reached, a value constrained by the total crystal content of typical Mamaku pumice. The amount of crystals in rhyolite-MELTS is determined by dividing the mass of the crystal phases by the mass of the system (excluding exsolved water). Separate isobaric cooling sequences are also calculated for the whole-rock compositions shown in Table 1, using the pressure retrieved from rhyolite-MELTS geobarometry of matrix glass (see Table 2). For each composition, sequences are calculated from initial conditions with water fractions of 0.1, 0.5, 1, 1.5, and 2 relative to saturation to allow for variability and uncertainty in the initial saturation conditions of the magma. Calculations are also allowed to continue until a threshold value of 10 wt% crystals is reached. The enthalpy values and mass calculated during the isobaric cooling sequences are used to calculate the total change in heat as a result of crystallization, ΔHtotal. 

### **6.3. Monte Carlo Modeling for Groundwater Flow** 

Uncertainty in magma storage pressures comes from uncertainty in the compositions of Mamaku glass compositions. In order to capture this, 300 compositions were generated from a Gaussian distribution of one of the Mamaku glass compositions from a prescribed standard deviation and mean, shown in Table 2. GL2306-G was chosen as the mean, and it has a corresponding standard deviation calculated from the multiple analyses used to calculate the mean composition. Each of these compositions was used as input for the rhyolite-MELTS geobarometer (Gualda & Ghiorso, 2014), 298 of which returned results, giving 298 storage pressures, of which the mean and standard deviation were calculated and used to generate a Gaussian distribution. The distribution of storage pressures was used to calculate a distribution of storage depths used in the Monte Carlo analysis of the Peclet number. 

A separate Monte Carlo analysis of the Peclet number  was conducted to account for uncertainties of the various parameters used in its calculation. Each of the parameters of  in Equation 8 have some natural variation, such as the distribution of the storage depth L, which was generated above. The ranges and distributions of all parameters used in the calculation of  are listed in Table 3. In the case of parameters for which uncertainties are well defined (e.g., storage pressure), we assumed a Gaussian distribution with known values for the standard deviation. In cases in which a distribution is not known, a uniform distribution was assumed in line with the Principle of Indifference (Keynes, 1921). The distributions for each parameter are assumed to be independent. For the Monte Carlo analysis, 10,000 random numbers were drawn from the 

14 of 17 

KELLY ET AL. 

**Geochemistry, Geophysics, Geosystems** 

10.1029/2021GC009794 

distribution of each parameter, the Peclet number was calculated for each set of random numbers, and then the resulting flux and heat loss timescale were calculated. 

The distribution of the Darcy velocity vz incorporated as input to the Monte Carlo Analysis of  was calculated by solving Equation 8 for vz 

and then using the following parameters from the present-day TVZ: c0  4.2  103J, 0  800  900 kg m3,   2.3 W m1 K1, L  7 km,  10   50 (Dempsey, 2012). The parameters listed are fixed and not chosen from a distribution, and they are assumed independent. The range of Darcy velocities is then 10 9 1 _vz_<sup></sup> 8 7.  10 to 4 8.  10 _m s_ , with an assumed uniform distribution. Darcy velocities and values of  prior to the eruption are difficult to constrain, since the speed of groundwater in the past cannot be calculated. The characteristics of the modern TVZ were used in modeling to give an idea of the minimum behavior, assuming there is not a large silicic magma body currently being stored in the crust that may impact Darcy velocity and . 

## **Conflict of Interest** 

The authors declare no conflicts of interest relevant to this study. 

## **Data Availability Statement** 

Supporting data sets and code are located in a repository on Open Science Framework (OSF). The repository can be found here: https://doi.org/10.17605/OSF.IO/MVGPW. 

##### **Acknowledgments** 

The authors would like to thank the MESSY group at Vanderbilt University for continuous discussion and input at various stages of the manuscript; David Furbish and Kristen Fauria for input at the latest stages; Vanderbilt's Earth and Environmental Sciences Honors Thesis Group 2019–2020; Lydia Harmon; Laura Blackstone; Mariah Jenkins, Betony Langenfeld, and Kate Pippenger; Sarah Smithies; Cole Burchiel; Frontiers Abroad for a great research experience that inspired this manuscript. This work was supported by Alberstadt-Reesman-Stearns Field Studies Funds from the Department of Earth and Environmental Sciences at Vanderbilt University. Comments by three reviewers and by journal editor Marie Edmonds are greatly appreciated. 

## **References** 

- Allan, A. S. R., Morgan, D. J., Wilson, C. J. N., & Millet, M. A. (2013). From mush to eruption in centuries: Assembly of the super-sized Oruanui magma body. _Contributions to Mineralogy and Petrology_ , _166_ , 143–164. https://doi.org/10.1007/s00410-013-0869-2 

- Andersen, C., & Weis, P. (2020). Heat transfer from convecting magma reservoirs to hydrothermal fluid flow systems constrained by coupled numerical modeling. _Geophysical Research Letters_ , _47_ (23). https://doi.org/10.1029/2020GL089463 

- Annen, C., Blundy, J. D., Leuthold, J., & Sparks, R. S. J. (2015). Construction and evolution of igneous bodies: Towards an integrated perspective of crustal magmatism. _Lithos_ , _230_ , 206–221. https://doi.org/10.1016/j.lithos.2015.05.008 

- Ardid, A., Archer, R., Bertrand, E., Sepulveda, F., Tarits, P., & Dempsey, D. (2021). Heat transfer through the Wairakei-Tauhara geothermal system quantified by multi-channel data modeling. _Geophysical Research Letters_ , _48_ (8), 1–10. https://doi.org/10.1029/2020GL092056 

- Asimow, P. D., & Ghiorso, M. S. (1998). Algorithmic modifications extending MELTS to calculate subsolidus phase relations. _American Mineralogist_ , _83_ (9–10), 1127–1132. https://doi.org/10.2138/am-1998-9-1022 

- Audétat, A., Miyajima, N., Wiesner, D., & Audinot, J.-N. (2021). Confirmation of slow Ti diffusion in quartz by diffusion couple experiments and evidence from natural samples. _Geology_ , _49_ (8), 963–967. https://doi.org/10.1130/g48785.1 

- Barker, S. J., Wilson, C. J. N., Illsley-Kemp, F., Leonard, G. S., Mestel, E. R. H., Mauriohooho, K., & Charlier, B. L. A. (2021). Taupō: An overview of New Zealand's youngest supervolcano. _New Zealand Journal of Geology and Geophysics_ , _64_ (2–3), 320–346. https://doi.org/ 10.1080/00288306.2020.1792515 

- Bégué, F., Deering, C. D., Gravley, D. M., Kennedy, B. M., Chambefort, I., Gualda, G. A. R., & Bachmann, O. (2014). Extraction, storage and eruption of multiple isolated magma batches in the paired Mamaku and Ohakuri eruption, Taupo Volcanic Zone, New Zealand. _Journal of Petrology_ , _55_ (8), 1653–1684. https://doi.org/10.1093/petrology/egu038 

- Bibby, H. M., Bennie, S. L., Stagpoole, V. M., & Caldwell, T. G. (1994). Resistivity structure of the Waimangu, Waiotapu, Waikite and Reporoa geothermal areas, New Zealand. _Geothermics_ , _23_ (5–6), 445–471. https://doi.org/10.1016/0375-6505(94)90013-2 

- Bibby, H. M., Caldwell, T. G., Davey, F. J., & Webb, T. H. (1995). Geophysical evidence on the structure of the Taupo Volcanic Zone and its hy- 

drothermal circulation. _Journal of Volcanology and Geothermal Research_ , _68_ (1–3), 29–58. https://doi.org/10.1016/0377-0273(95)00007-H Boden, D. R. (2016). Geology and heat architecture of the Earth's interior. In _Geologic fundamentals of geothermal energy_ (pp. 53–70). CRC Press. https://doi.org/10.1201/9781315371436-4 

- Bredehoeft, J. D., & Papadapulos, I. S. (1965). Rates of vertical groundwater movement estimated from the Earth's thermal profile. _Water Resources Research_ , _1_ (2), 325–328. https://doi.org/10.1029/WR001i002p00325 

- Bryan, S. E., Peate, I. U., Peate, D. W., Self, S., Jerram, D. A., Mawby, M. R., et al. (2010). The largest volcanic eruptions on Earth. _Earth-Science Reviews_ , _102_ (3–4), 207–229. https://doi.org/10.1016/j.earscirev.2010.07.001 

- Carslaw, H. S., & Jaeger, J. C. (1984). _Conduction of heat in solids_ . Oxford University Press. 

- Cashman, K. V., & Giordano, G. (2014). Calderas and magma reservoirs. _Journal of Volcanology and Geothermal Research_ , _288_ , 28–45. https://doi.org/10.1016/j.jvolgeores.2014.09.007 

- Cathles, L. M., Erendi, A. H. J., & Barrie, T. (1997). How long can a hydrothermal system be sustained by a single intrusive event. _Economic Geology_ , _92_ (7–8), 766–771. https://doi.org/10.2113/gsecongeo.92.7-8.766 

- Charlier, B. L. A., Wilson, C. J. N., & Davidson, J. P. (2008). Rapid open-system assembly of a large silicic magma body: Time-resolved evidence from cored plagioclase crystals in the Oruanui eruption deposits, New Zealand. _Contributions to Mineralogy and Petrology_ , _156_ (6), 799–813. https://doi.org/10.1007/s00410-008-0316-y 

15 of 17 

KELLY ET AL. 

**Geochemistry, Geophysics, Geosystems** 

10.1029/2021GC009794 

Cooper, G. F., Wilson, C. J. N., Millet, M. A., Baker, J. A., & Smith, E. G. C. (2012). Systematic tapping of independent magma chambers during the 1Ma Kidnappers supereruption. _Earth and Planetary Science Letters_ , _313_ (314), 23–33. https://doi.org/10.1016/j.epsl.2011.11.006 Cooper, K. M., & Kent, A. J. R. (2014). Rapid remobilization of magmatic crystals kept in cold storage. _Nature_ , _506_ , 480–483. https://doi. org/10.1038/nature12991 

Dempsey, D. E. (2012). _Numerical models of tectonism and geothermal circulation with application to the Taupo Volcanic Zone, New Zealand_ (Doctoral dissertation). University of Auckland. https://doi.org/10.5281/zenodo.4793934 

Dempsey, D. E., Gravley, D. M., & Rowland, J. (2018). _Evolution and eruption of geothermally cooled magma bodies_ . Paper presented at 40th New Zealand Geothermal Workshop, New Zealand Geothermal Workshop. https://doi.org/10.5281/zenodo.4794066 

Dempsey, D. E., Simmons, S. F., Archer, R. A., & Rowland, J. V. (2012). Delineation of catchment zones of geothermal systems in largescale rifted settings. _Journal of Geophysical Research_ , _117_ (B10). https://doi.org/10.1029/2012JB009515 

de Silva, S. L., & Gosnold, W. D. (2007). Episodic construction of batholiths: Insights from the spatiotemporal development of an ignimbrite flare-up. _Journal of Volcanology and Geothermal Research_ , _167_ (1–4), 320–335. https://doi.org/10.1016/j.jvolgeores.2007.07.015 

Druitt, T. H., Costa, F., Deloule, E., Dungan, M., & Scaillet, B. (2012). Decadal to monthly timescales of magma transfer and reservoir growth at a caldera volcano. _Nature_ , _482_ (7383), 77–80. https://doi.org/10.1038/nature10706 

Elston, W. E. (1984). Mid-tertiary ash flow tuff cauldrons, southwestern New Mexico. _Journal of Geophysical Research_ , _89_ (B10), 8733–8750. https://doi.org/10.1029/JB089iB10p08733 

Eppelbaum, L., Kutasov, I., & Pilchin, A. (2014). Thermal properties of rocks and density of fluids. In _Applied geothermics_ (pp. 99–149). Springer-Verlag. https://doi.org/10.1007/978-3-642-34023-9_2 

Fournier, R. O. (1999). Hydrothermal processes related to movement of fluid from plastic into brittle rock in the magmatic-epithermal environment. _Economic Geology_ , _94_ (8), 1193–1211. https://doi.org/10.2113/gsecongeo.94.8.1193 

Ghiorso, M. S., & Sack, R. O. (1995). Chemical mass transfer in magmatic processes IV. A revised and internally consistent thermodynamic model for the interpolation and extrapolation of liquid-solid equilibria in magmatic systems at elevated temperatures and pressures. _Contributions to Mineralogy and Petrology_ , _119_ (2), 197–212. https://doi.org/10.1007/BF00307281 

Giggenbach, W. F. (1995). Variations in the chemical and isotopic composition of fluids discharged from the Taupo Volcanic Zone, New Zealand. _Journal of Volcanology and Geothermal Research_ , _68_ (1–3), 89–116. https://doi.org/10.1016/0377-0273(95)00009-J 

Gravley, D. M., Deering, C. D., Leonard, G. S., & Rowland, J. V. (2016). Ignimbrite flare-ups and their drivers: A New Zealand perspective. _Earth-Science Reviews_ , _162_ , 65–82. https://doi.org/10.1016/j.earscirev.2016.09.007 

Gravley, D. M., Wilson, C. J. N., Leonard, G. S., & Cole, J. W. (2007). Double trouble: Paired ignimbrite eruptions and collateral subsidence in the Taupo Volcanic Zone, New Zealand. _GSA Bulletin_ , _119_ (1–2), 18–30. https://doi.org/10.1130/B25924.1 Gualda, G. A. R., & Ghiorso, M. S. (2013a). Low-pressure origin of high-silica rhyolites and granites. _Journal of Geology_ , _121_ (5), 537–545. https://doi.org/10.1086/671395 

Gualda, G. A. R., & Ghiorso, M. S. (2013b). The Bishop Tuff giant magma body: An alternative to the Standard Model. _Contributions to Mineralogy and Petrology_ , _166_ (3), 755–775. https://doi.org/10.1007/s00410-013-0901-6 

Gualda, G. A. R., & Ghiorso, M. S. (2014). Phase-equilibrium geobarometers for silicic rocks based on rhyolite-MELTS. Part 1: Principles, procedures, and evaluation of the method. _Contributions to Mineralogy and Petrology_ , _168_ , 1033. https://doi.org/10.1007/ s00410-014-1033-3 

Gualda, G. A. R., Ghiorso, M. S., Lemons, R. V., & Carley, T. L. (2012). Rhyolite-MELTS: A modified calibration of MELTS optimized for silica-rich, fluid-bearing magmatic systems. _Journal of Petrology_ , _53_ (5), 875–890. https://doi.org/10.1093/petrology/egr080 Gualda, G. A. R., Gravley, D. M., Connor, M., Hollmann, B., Pamukcu, A. S., Bégué, F., et al. (2018). Climbing the crustal ladder: Magma storage-depth evolution during a volcanic flare-up. _Science Advances_ , _4_ (10), eaap7567. https://doi.org/10.1126/sciadv.aap7567 

Gualda, G. A. R., Gravley, D. M., Deering, C. D., & Ghiorso, M. S. (2019). Magma extraction pressures and the architecture of volcanic plumbing systems. _Earth and Planetary Science Letters_ , _522_ , 118–124. https://doi.org/10.1016/j.epsl.2019.06.020 

Gualda, G. A. R., & Pamukçu, A. S. (2020). New Ti-in-quartz diffusivites reconcile natural Ti zoning with time scales and temperature of upper crustal magma reservoirs: Comment. _Geology_ , _48_ (12), e513. https://doi.org/10.1130/G48242C.1 

Gualda, G. A. R., Pamukcu, A. S., Ghiorso, M. S., Anderson, A. T., Sutton, S. R., & Rivers, M. L. (2012). Timescales of Quartz crystallization and the longevity of the Bishop Giant Magma body. _PLoS One_ , _7_ (5), e37492. https://doi.org/10.1371/journal.pone.0037492 

Gualda, G. A. R., & Sutton, S. R. (2016). The year leading to a supereruption. _PLoS One_ , _11_ (7), e0159200. https://doi.org/10.1371/journal. pone.0159200 Hardee, H. C. (1981). Convective heat extraction from molten magma. _Journal of Volcanology and Geothermal Research_ , _10_ (1–3), 175–193. https://doi.org/10.1016/0377-0273(81)90061-5 

Hardee, H. C., & Dunn, J. C. (1981). Convective heat transfer in magmas near the liquidus. _Journal of Volcanology and Geothermal Research_ , _10_ (1–3), 195–207. https://doi.org/10.1016/0377-0273(81)90062-7 

Hardee, H. C., & Larson, D. W. (1977). The extraction of heat from magmas based on heat transfer mechanisms. _Journal of Volcanology and Geothermal Research_ , _2_ (2), 113–144. https://doi.org/10.1016/0377-0273(77)90006-3 Hayba, D. O., & Ingebritsen, S. E. (1997). Multiphase groundwater flow near cooling plutons. _Journal of Geophysical Research_ , _102_ (B6), 12235–12252. https://doi.org/10.1029/97jb00552 Hedenquist, J. W. (1990). The thermal and geochemical structure of the Broadlands-Ohaaki geothermal system, New Zealand. _Geothermics_ , _19_ (2), 151–185. https://doi.org/10.1016/0375-6505(90)90014-3 

Hochstein, M. P. (1995). Crustal heat transfer in the Taupo Volcanic Zone (New Zealand): Comparison with other volcanic arcs and explanatory heat source models. _Journal of Volcanology and Geothermal Research_ , _68_ (1–3), 117–151. https://doi.org/10.1016/0377-0273(95)00010-R Houghton, B. F., Wilson, C. J. N., McWilliams, M. O., Lanphere, M. A., Weaver, S. D., Briggs, R. M., & Pringle, M. S. (1995). Chronology and dynamics of a large silicic magmatic system: Central Taupo Volcanic Zone, New Zealand. _Geology_ , _23_ (1), 13–16. https://doi. org/10.1130/0091-7613(1995)023<0013:CADOAL>2.3.CO;2 

Jollands, M. C., Bloch, E., & Müntener, O. (2020a). New Ti-in-quartz diffusivities reconcile natural Ti zoning with time scales and temperatures of upper crustal magma reservoirs. _Geology_ , _48_ (7), 654–657. https://doi.org/10.1130/G47238.1 

Jollands, M. C., Bloch, E., & Müntener, O. (2020b). New Ti-in-quartz diffusivities reconcile natural Ti zoning with time scales and temperatures of upper crustal magma reservoirs: Reply. _Geology_ , _48_ (12), e514. https://doi.org/10.1130/G48470Y.1 

Keynes, J. M. (1921). The principle of indifference. In _A treatise on probability_ (pp. 44–70). MacMillan & Co. 

Kissling, W. M., & Weir, G. J. (2005). The spatial distribution of the geothermal fields in the Taupo Volcanic Zone, New Zealand. _Journal of Volcanology and Geothermal Research_ , _145_ (1–2), 136–150. https://doi.org/10.1016/j.jvolgeores.2005.01.006 

Lipman, P. W. (1984). The roots of ash flow calderas in western North America: Windows into the tops of granitic batholiths. _Journal of Geophysical Research_ , _89_ (B10), 8801–8841. https://doi.org/10.1029/JB089iB10p08801 

16 of 17 

KELLY ET AL. 

**Geochemistry, Geophysics, Geosystems** 

10.1029/2021GC009794 

Milner, D. M., Cole, J. W., & Wood, C. P. (2002). Asymmetric, multiple-block collapse at Rotorua Caldera, Taupo Volcanic Zone, New Zealand. _Bulletin of Volcanology_ , _64_ (2), 134–149. https://doi.org/10.1007/s00445-001-0191-0 

Mittal, T., & Richards, M. A. (2019). Volatile degassing from magma chambers as a control on volcanic eruptions. _Journal of Geophysical Research: Solid Earth_ , _124_ (8), 7869–7901. https://doi.org/10.1029/2018JB016983 

Norton, D., & Taylor, H. P. (1979). Quantitative simulation of the hydrothermal systems of crystallizing magmas on the basis of transport theory and oxygen isotope data: An analysis of the skaergaard intrusion. _Journal of Petrology_ , _20_ (3), 421–486. https://doi.org/10.1093/ petrology/20.3.421 

O'Sullivan, M. J., Yeh, A., & Mannington, W. I. (2009). A history of numerical modelling of the Wairakei geothermal field. _Geothermics_ , _38_ (1), 155–168. https://doi.org/10.1016/j.geothermics.2008.12.001 

Pamukcu, A. S., Gualda, G. A. R., Bégué, F., & Gravley, D. M. (2015). Melt inclusion shapes: Timekeepers of short-lived giant magma bodies. _Geology_ , _43_ (11), 947–950. https://doi.org/10.1130/G37021.1 

Pamukcu, A. S., Gualda, G. A. R., Ghiorso, M. S., Miller, C. F., & McCracken, R. G. (2015). Phase-equilibrium geobarometers for silicic rocks based on rhyolite-MELTS—Part 3: Application to the Peach Spring Tuff (Arizona–California–Nevada, USA). _Contributions to Mineralogy and Petrology_ , _169_ (3), 33. https://doi.org/10.1007/s00410-015-1122-y 

Perrot, P. (1998). _A to Z of thermodynamics_ . Oxford University Press. 

Pitcher, B. W., Gualda, G. A. R., & Hasegawa, T. (2020). Repetitive duality of rhyolite compositions, timescales, and storage and extraction conditions for pleistocene caldera-forming eruptions, Hokkaido, Japan. _Journal of Petrology_ , _62_ (2), egaa106. https://doi.org/10.1093/ petrology/egaa106 

Rowland, J. V., & Sibson, R. H. (2004). Structural controls on hydrothermal flow in a segmented rift system, Taupo Volcanic Zone, New Zealand. _Geofluids_ , _4_ , 259–283. https://doi.org/10.1111/j.1468-8123.2004.00091.x 

Rowland, J. V., & Simmons, S. F. (2012). Hydrologic, magmatic, and tectonic controls on hydrothermal flow, Taupo Volcanic Zone, New Zealand: Implications for the formation of epithermal vein deposits. _Economic Geology_ , _107_ (3), 427–457. https://doi.org/10.2113/ econgeo.107.3.427 

Scholz, C. H. (1988). The brittle-plastic transition and the depth of seismic faulting. _Geologische Rundschau_ , _77_ (1), 319–328. https://doi. org/10.1007/BF01848693 

Shamloo, H. I., & Till, C. B. (2019). Decadal transition from quiescence to supereruption: Petrologic investigation of the Lava Creek Tuff, Yellowstone Caldera, WY. _Contributions to Mineralogy and Petrology_ , _174_ (4), 32. https://doi.org/10.1007/s00410-019-1570-x 

Weir, G. J. (2009). A mathematical model of rainfall-controlled geothermal fields. _Transport in Porous Media_ , _77_ (2), 323–334. https://doi. org/10.1007/s11242-008-9312-y 

Wilson, C. J. N., Houghton, B. F., McWilliams, M. O., Lanphere, M. A., Weaver, S. D., & Briggs, R. M. (1995). Volcanic and structural evolution of Taupo Volcanic Zone, New Zealand: A review. _Journal of Volcanology and Geothermal Research_ , _68_ (1–3), 1–28. https://doi. org/10.1016/0377-0273(95)00006-G 

Wilson, C. J. N., & Rowland, J. V. (2016). The volcanic, magmatic and tectonic setting of the Taupo Volcanic Zone, New Zealand, reviewed from a geothermal perspective. _Geothermics_ , _59_ (Part B), 168–187. https://doi.org/10.1016/j.geothermics.2015.06.013 

17 of 17 

KELLY ET AL.