Catastrophic caldera‐forming volcanoes, often referred to as “supervolcanoes,” are one of the greatest natural hazards on Earth. Understanding how these explosive eruptions are triggered is critical for forecasting volcanic activity at the world's largest volcanoes and mitigating their hazards. Using sophisticated numerical models, the affect of plate tectonic stress is investigated to determine its role in triggering large caldera‐forming eruptions. The numerical experiments indicate that calderas located in extensional tectonic settings erupt more readily than calderas located in compressional settings or in tectonically neutral settings. Additionally, calderas in extensional settings fail on shorter timescales than calderas located in compressional or tectonically neutral settings. The most important outcome of this work is that our numerical models show for the first time that the rock surrounding a large magma reservoir is only stable on timescales of centuries to thousands of years when new magma is actively being injected into the magma reservoir. This finding provides important constraints on the amount of time necessary to recharge and erupt a large, supervolcano size reservoir from the first indication of magmatic activity.

We utilize 3‐D temperature‐dependent viscoelastic finite element models to investigate the mechanical response of the host rock supporting large caldera‐size magma reservoirs (volumes >10 2 km 3 ) to local tectonic stresses. The mechanical stability of the host rock is used to determine the maximum predicted repose intervals and magma flux rates that systems may experience before successive eruption is triggered. Numerical results indicate that regional extension decreases the stability of the roof rock overlying a magma reservoir, thereby promoting early‐onset caldera collapse. Alternatively, moderate amounts of compression (≤10 mm/year) on relatively short timescales (<10 4 years) increases roof rock stability. In addition to quantifying the affect of tectonic stresses on reservoir stability, our models indicate that the process of rejuvenation and mechanical failure is likely to take place over short time periods of hundreds to thousands of years. These findings support the short preeruption melt accumulation timescales indicated by U series disequilibrium studies.

1 Introduction The classic paradigm in volcanology is that eruption commences when the pressure within a magma reservoir exceeds the confining strength of the surrounding host rock (Sparks, 2003). However, the link between overpressure and eruption becomes tenuous as reservoirs reach sizes >103 km3. At such scales, the heat flux from what are presumably long‐lived magmatic systems greatly impacts the rheological properties of the surrounding rock, limiting the potential for magma injection or other changes to elevate internal overpressure and thus buffering large magma reservoirs from failure (Gregg et al., 2012; Jellinek & DePaolo, 2003). In the absence of an internal overpressure‐driven mechanism, faulting and tectonic forcing has been suggested as a catalyst for rupturing large, shallow magma bodies to trigger caldera‐forming eruptions (Allan et al., 2012; Gregg et al., 2012, 2015). In this scenario the mechanical stability of the host rock is critical, and local tectonic stress may have a profound affect on the corresponding stability of large magma reservoirs. Previous investigations suggest that regional extension may promote the growth of large silicic magma reservoirs within the upper crust by accommodating the accumulation of large volumes of melt (Jellinek & DePaolo, 2003). Indeed, statistical analysis of 108 large, Quaternary, silicic calderas indicates that approximately 70% of calderas are located in extensional or transtensional tectonic regimes (Figure S1 in the supporting information; Hughes & Mahood, 2008, 2011). Additionally, recent investigations from the Taupo Volcanic Zone (TVZ) suggest tectonic controls on eruptive products and processes (Spinks et al., 2005) and on the mechanics of triggering the Oruanui eruption at Taupo Caldera (Allan et al., 2012). Furthermore, numerical results (e.g., Currenti & Williams, 2014; Geyer & Marti, 2009; Gray & Monaghan, 2004; Jellinek & DePaolo, 2003) and analog models (Acocella, 2007, 2010) have indicated that extensional tectonic regimes promote fault development and failure. Given the well‐established link between tectonic regime and propensity for developing large caldera systems, it is vital to quantitatively assess the timescales at which a magma reservoir becomes unstable under tectonically active conditions. In this study we utilize three‐dimensional, temperature‐dependent, viscoelastic finite element models (FEMs), which incorporate a Mohr‐Coulomb failure criterion, to investigate the role of regional tectonic stress as a trigger for caldera‐forming eruptions (Gregg et al., 2012; Grosfils et al., 2015; see Supplemental Methods section in the supporting information for details). Numerical experiments are conducted to systematically investigate magma reservoir stability for a global range of far‐field stresses and magmatic flux rates. Maximum repose timescales are calculated to determine the impact of tectonic stress and flux rates on magma chamber stability at long timescales. We apply our model to Taupo Caldera, New Zealand, as a case study to test model‐predicted repose times against those observed. Results indicate that tectonic regime plays an integral role, impacting eruption repose times and both the longevity and stability of large upper crustal magma systems. This study provides the first quantitative approach to address the mechanical relationship between local tectonic stress and the timescales of magma reservoir stability.

2 Tectonic Modulation of Caldera Eruptions Building upon the two‐dimensional approach of Gregg et al. (2012), five 3‐D, thermomechanical numerical experiments are conducted to investigate the relationship between tectonic stress and reservoir stability. (See Supplementary Methods section in the supporting information for model details.) Three numerical experiments investigate magma reservoir stability for a neutral tectonic regime, an extensive regime, and a compressive regime at a range of magma flux rates. For comparison purposes, two additional numerical experiments determine reservoir stability in an extensional and compressional setting when there is no influx of new magma. The stability of the host rock is determined by predictions of Mohr‐Coulomb failure. When the FEM predicts throughgoing fault development within the roof overlying the magma reservoir (i.e., Mohr‐Coulomb failure is calculated from the surface of the model space to the reservoir boundary), the magma reservoir is no longer considered stable. The timing required for throughgoing fault development indicates the maximum repose time for a given set of model conditions. The imposed tectonic stresses and magma flux rates are maintained as constant throughout the numerical experiment to provide an end‐member approximation. We calculate repose times of magma reservoirs located in tectonically neutral regimes, expanding at a range of flux rates, to establish a standard to compare with reservoirs located in tectonically active regimes (Figure 1). Applying a constant rate of applied extensional stress to the host rock surrounding a reservoir that is expanding yields repose times that are shorter than for a reservoir in a neutral tectonic setting. In addition, the maximum flux rates that can be accommodated before an eruption is triggered are lower for a reservoir in an extensional tectonic setting in every case. Conversely, the maximum flux rates that can be accommodated are higher for reservoirs subject to a constant rate of applied compressive stress. This relationship occurs in each initial reservoir volume tested for this study, and we observe that maximum flux rates and reservoir growth are correlated with initial reservoir volume. Figure 2 illustrates this relationship by comparing reservoir volumes, ranging from 102 to 105 km3, which are subject to no external tectonic stress (Figure 2). This finding suggests similar mechanical failure triggers among a variety of caldera sizes. Figure 1 Open in figure viewer PowerPoint Maximum model‐predicted flux rates for reservoir volumes of 103 and 104 km3. (a and c) The maximum model‐predicted flux rates and (b and d) the maximum model‐predicted reservoir volume change as a percent change from the initial volume. The solid black line indicates the maximum possible flux rate before failure is initiated in a system experiencing no external stress. Using this line as a reference, it is apparent that systems experiencing compressive stress (red) can withstand higher fluxes than those not exposed to external stress. Systems experiencing extensive stress (blue) withstand lower flux rates than those in compressive or neutral tectonic regimes. The green line represents the case of Taupo Caldera, which is located in an extensive regime characterized by 7.2‐mm/year extension. Reservoirs fail with no added volumetric flux at the time data collection was terminated. Note the change in scale for the case of the 104 km3 reservoir. Figure 2 Open in figure viewer PowerPoint Maximum flux rates required to trigger failure, plotted for a variety initial reservoir volumes. The plotted values are for systems subject to no external tectonic stress. Maximum flux rates scale with the initial reservoir volumes. For the purpose of understanding the dependence of reservoir stability upon external stress alone, we assess reservoirs that are not expanding but are subject only to a constant rate of applied external stress. A reservoir that is expanding due to either magma flux or volatile exsolution will always fail before a system, which is dormant. Consequently, our nonexpanding case studies constrain the maximum repose times for calderas located in tectonically active regimes, assuming that the reservoirs are thermally primed for eruption. Results indicate that repose times for calderas located in compressive regimes are longer than those of calderas located in extensional regimes characterized by the same magnitude of stress. However, repose times for calderas located in compressive regimes are still shorter than those of calderas located in tectonically neutral regimes. This finding suggests that while regional compression may temporarily strengthen the host rock and buffer it from failure, on longer timescales compressional stress weakens the system and may promote a caldera‐forming eruption earlier than is forecasted for systems without any external tectonic stress. Similarly, this finding yields support to a recently proposed hypothesis that following periods of repose and rapid inflation, stresses within the crust are not alleviated but are retained. In this way, stresses build until energy is released in a series of faulting events, and eruptions are triggered when bulk failure of the host rock is achieved (Kilburn et al., 2017). It is important to note that reservoirs that are not being rejuvenated through recharge will likely cool and crystallize before the critical timescale is reached; however, these results illustrate the importance of considering external stress when evaluating the stability of a large magma reservoir system. Previous thermomechanical models have demonstrated that, for the largest caldera systems (reservoir volumes >103 km3), eruptions are likely triggered when flexure of the roof overlying a magma reservoir generates downward propagating faults that facilitate rapid magma evacuation and caldera collapse (Gregg et al., 2012, 2015). When coupled with the added effect of external tectonic stress, the flexing roof becomes even weaker and more prone to early onset collapse. Extending the roof overlying an expanding magma reservoir causes stress to concentrate in the host rock above the apex of the reservoir. Consequently, throughgoing fault development in extensional regimes is observed in the overlying roof rock where vertical surface deformation is extreme and the roof rock is thin (Figure 3). The rapid generation of ring faults in the roof causes roof destabilization and ultimately triggers eruption. Alternatively, applying regional compression temporarily strengthens and stabilizes the roof overlying a reservoir by concentrating stress at the distal edges above a reservoir where surface deformation is less intense and the crust is relatively thick. This allows reservoirs to accommodate higher fluxes of material and greater reservoir volume changes before throughgoing faults develop and subsequently trigger a caldera‐forming eruption. Figure 3 Open in figure viewer PowerPoint The 2‐D slice through the 3‐D model outputs plot total stress only where failure occurs for a 1 km3 reservoir at 0, 5,000 and 10,000 years. Caldera collapse is defined as throughgoing Mohr‐Coulomb failure and is observed at 5,000 years given an extension rate of 7.2 and 20 mm/year. Throughgoing Mohr‐Coulomb failure is not observed at 5,000 years given the case of a convergence rate of 20 mm/year but is observed later at 10,000 years. The numerical models assume a pristine, homogenous host rock, unlikely in a natural setting since structural flaws and other shallow crustal heterogeneities will exist in volcanically and tectonically active regions (Geyer & Marti, 2009; Hughes & Mahood, 2011). Preexisting weaknesses will act to reduce repose intervals and maximum flux rates. Furthermore, the models do not consider periodicity of magmatic recharge and inflation or dormancy and subsidence as is observed at most calderas. Likewise, small eruptions occurring between caldera‐forming events may affect the distribution and magnitude of stress in the crust and are not considered in this study. Kilburn et al. (2017) recently proposed that stresses building in the crust during such events may not be alleviated during dormancy and may play a significant role in triggering future eruptions. As such, our calculated repose intervals and maximum magma flux rates provide end‐member approximations.

3 Case Study: Forecasting Repose Intervals for Taupo Caldera Located in a regionally extensional back‐arc spreading center, the TVZ is one of the most productive volcanic regions on Earth (Wilson et al., 1995). The largest of the calderas located within the TVZ, Taupo Caldera, is located in an area of the spreading center characterized by 7.2 mm/year of pure extension with no dextral shear stress (Spinks et al., 2005; Figure 4). Since the establishment of the major structural features defining the young TVZ 350 ka, two caldera collapse events have been documented at Taupo Caldera—the caldera‐forming Oruanui eruption (26.5 ka, 530 km3) and the Taupo eruption (1.8 ka, 35 km3), with 27 minor eruptions occurring in interim (Barker et al., 2015; Houghton et al., 1995; Wilson & Charlier, 2009; Wilson & Rowland, 2016). The minor eruptions occurring between these two events collectively account for 28.25 km3 of extruded materials, and thus, only 2.4% of the total bulk volume erupted from Taupo Caldera in the past thirty thousand years (Barker et al., 2015; Wilson, 1993). While the smaller eruptions of Taupo Caldera may temporarily impact the host rock stress state, we assume given their localized nature and small magnitudes that they are of minimal long‐term mechanical importance (a hypothesis that would benefit from future study), and they are not considered for this study of large caldera‐forming eruptions. As such, we define the caldera repose interval as the amount of time between two caldera‐forming eruptions (the repose interval between the two collapse events at Taupo Caldera is 24.7 ky). Although we test a range of initial reservoir volumes in our parameterized experiments, for the Taupo application we assume 10% reservoir evacuation (e.g., Gregg et al., 2012) to estimate initial reservoir volumes for the Oruanui and Taupo eruptions. Figure 4 Open in figure viewer PowerPoint 3 km3 reservoir with the overlain topography of Taupo Caldera. The reservoir (outlined in red) is treated as a pressurized cavity embedded within a temperature‐dependent viscoelastic medium. Boundary conditions on two of the model‐bounding surfaces are prescribed with velocities to simulate the affect of external tectonic stress. Arrows indicate the direction of spreading within the Taupo Volcanic Zone and consequently that used in the FEM. The structural caldera as defined by Wilson and Charlier ( 2009 3 km3 reservoir size. Model set‐up showing a 10kmreservoir with the overlain topography of Taupo Caldera. The reservoir (outlined in red) is treated as a pressurized cavity embedded within a temperature‐dependent viscoelastic medium. Boundary conditions on two of the model‐bounding surfaces are prescribed with velocities to simulate the affect of external tectonic stress. Arrows indicate the direction of spreading within the Taupo Volcanic Zone and consequently that used in the FEM. The structural caldera as defined by Wilson and Charlier () is superimposed on the model's surface (black dashed line), showing a good fit with the 10kmreservoir size. Figure 5 provides an estimate of the maximum repose interval and reservoir volume combinations given a constant stress regime and magma injection rate. Magmatic flux rates of 10−4, 10−3, 10−2, and 10−1 km3/year are tested to constrain repose intervals and reservoir volume combinations promoting failure; these modeled flux rates coincide with the previously suggested values of Wilson et al. (2006), Annen (2009), and Gelman et al. (2013) of 1.3E‐2, 1E‐2, and 5E‐3–8E‐3 km3/year, respectively. At these flux rates, model‐predicted repose stability periods at Taupo Caldera are far shorter than observed, indicating either that flux rates are far lower or that constant recharge is unlikely. The latter interpretation is intuitive since it is unlikely that a magma reservoir will grow at a constant rate for the duration of its repose period, and indeed, variable flux has been proposed by several studies of Taupo Caldera (Allan et al., 2013, 2017; Barker et al., 2015). Our model results indicate that intermittent periods of high flux on relatively short timescales, separated by longer periods of stasis, produce the repose periods observed at Taupo Caldera. Figure 5 Open in figure viewer PowerPoint Predictions of repose interval for reservoirs ranging in size from 102 to 104 km3 given applied flux rates ranging in magnitude from 10−4 to 10−1 km3/year. Solid lines indicate systems, which are not exposed to any external stress, while dashed lines represent the tectonic regime characteristic of the Taupo caldera (7.2‐mm/year extension). Estimates of the pre‐Oruanui and the pre‐Taupo eruption reservoir volumes are indicated by the vertical dashed black lines and represent the size of the reservoir if the erupted volumes were 10% of the initial reservoir volumes. Model predictions of shorter‐than‐observed repose intervals indicate that magma reservoirs will not remain mechanically stable during prolonged periods of increased magma flux. This finding is consistent with recent paleouplift and resurgence studies suggesting surface deformation often occurs episodically during periods of quiescence rather than as one continuous event preceding a large eruption, suggesting intermittent periods of magmatic fluxing (Hickey et al., 2013; Mucek et al., 2017; Perkins et al., 2016). This mechanical result is supported by recent uranium series disequilibria studies indicating that magma reservoirs are not maintained in an easily eruptible state (>750 °C and below the critical crystallinity of 40–50%) for long intervals by periodic recharge of hotter materials but rather are remobilized from cold storage via a recharge event and near‐simultaneous eruption (Burgisser & Bergantz, 2011; Cooper & Kent, 2014). Furthermore, U‐Th disequilibrium model‐age observations from Oruanui ignimbrite samples, along with element diffusion modeling, suggest a transient thermal behavior of the Oruanui melt body. Specifically, results indicate that while the Oruanui mush source developed over tens of thousands of years, the eruptible reservoir formed in ≤3,000 years (Allan et al., 2013; Rubin et al., 2017; Wilson & Charlier, 2009). Our models utilize a pressurized cavity to simulate an idealized fluid and, as such, do not consider magmatic conditions such as compressibility explicitly. Furthermore, we simulate a thermally primed, eruptible reservoir. For this reason, our model results indicate the maximum repose times for a thermally primed system experiencing constant magmatic flux and tectonic conditions. The relatively short geochemical estimates of magma accumulation timescales are similar to those predicted for repose timescales in the FEM when conditions preceding eruption are constant. Therefore, our numerical results suggest from a mechanical perspective that the magma reservoirs feeding caldera systems remain stable for only short periods of time when experiencing high magmatic flux rates.

4 Discussion and Conclusions Three‐dimensional thermomechanical FEMs indicate that external tectonic stress plays a critical role in the timing and triggering of caldera‐forming eruptions. In particular, extensional tectonic stress regimes promote failure in the host rock leading to destabilization and eruption. Systems located in extensive regimes have shorter repose intervals and fail with lower magmatic fluxes than those located in compressional regimes. On short timescales, compressional tectonic stress appears to temporarily stabilize the roof rock surrounding a reservoir; however, on longer timescales compressional systems have shorter repose periods than systems in a neutral tectonic setting. This phenomenon is directly tied to the way in which strain is accumulated in the roof overlying an expanding reservoir. Model outputs are used to infer maximum repose times for given magmatic flux rates and suggest that magma reservoirs are stable only for hundreds to thousands of years while undergoing a magmatic fluxing event. These findings are consistent with melt accumulation timescales determined by U series disequilibrium studies and provide a physics‐grounded mechanical framework that directly supports rapid melt remobilization from cold storage (e.g., Cooper & Kent, 2014). Furthermore, model‐predicted repose timescales for the Oruanui eruption are consistent with the magma accumulation timescale suggested by Allan et al. (2013) and Rubin et al. (2017) for eruptions occurring within the TVZ. Findings of this study thus have broader implications for assessing eruption potential for calderas actively exhibiting signs of unrest.

Acknowledgments This work was supported by a University of Illinois startup fund (Cabaniss and Gregg), by NSF (OCE1634995, Cabaniss and Gregg) and by NASA (NNX12AO49G, Grosfils). We gratefully acknowledge helpful insightful reviews by G. De Natale, C. Wauthier, T. Masterlark, and an anonymous reviewer, which greatly improved this manuscript. We are grateful for helpful comments by and discussions with S. de Silva, C. Wilson, L. Karlstrom, D. Burns, D. Geist, B. Chadwick, S. Nooner, S. Marshak, Y. Zhan, J. Albright, R. Goldman, and the UIUC Geodynamics Group. Data in support of this manuscript are published in the PANGAEA data repository and available online (https://doi.pangaea.de/10.1594/PANGAEA.885992).

Supporting Information Filename Description grl57310-sup-0001-2018GL077393-S01.pdfPDF document, 3.2 MB Supporting Information S1 Please note: The publisher is not responsible for the content or functionality of any supporting information supplied by the authors. Any queries (other than missing content) should be directed to the corresponding author for the article.