Stretching, Shaking, Inflating: Volcanic-Tectonic Interactions at a Rifting Silicic Caldera

1 Introduction

Silicic (dacite to rhyolite compositions in this context) caldera volcanoes have demonstrated their capability to produce the largest and most catastrophic explosive volcanic eruptions on Earth. There is thus a critical need to constrain the various processes that lead to unrest or eruptions within these systems. To better understand the dynamics of silicic volcanoes, they are often examined through geochemical and/or physical volcanological lenses to address the mechanisms of growth of their magma bodies, magma ascent dynamics, and the physical processes controlling eruptions and their widespread tephra deposits (see reviews from Cashman and Sparks, 2013; Wilson et al., 2021). With perhaps the exception of the 1991 Pinatubo eruption (Philippines) (Harlow et al., 1991; Jones et al., 2001), the events leading up to caldera-forming silicic eruptions have not been monitored in detail during modern (instrumental) times, and signs of unrest vary widely among different active calderas (Acocella et al., 2015). As such, there are still large uncertainties over the deformation signals that would be observed at the surface prior to and during caldera-forming silicic eruption. However, deformation associated with past eruptive and non-eruptive unrest events is preserved to varied degrees in both the historical and geological records and can provide useful insights into the processes that operate at caldera volcanoes over a range of timescales.

Most silicic caldera volcanoes occur in areas of local extensional or transtensional stress (Hughes and Mahood, 2011), and episodes of volcanic and earthquake activity may be triggered and modulated by the various interactions between the magmatic (basaltic and felsic) and tectonic systems (Rowland et al., 2010). Indeed, the presence of faults in the crust and variable static and/or dynamic stress changes in extensional settings has been shown to impact magma migration and the eruptive behavior of a volcano (e.g., Buck et al., 2006; Baer et al., 2008; Menand and Daniels, 2010; Biggs et al., 2013; Le Corvec et al., 2013; Acocella et al., 2015; Gómez-Vasconcelos et al., 2020a). For example, the simultaneous and sequential tapping of discrete magma bodies during the 2.1 Ma Huckleberry Ridge eruption (Yellowstone) is interpreted to support a possible external tectonic trigger to this eruptive episode (Myers et al., 2016). Also, the 1350 CE Mono-Inyo eruptions in California were associated with silicic dike injection and shallow normal fault slip (Mastin and Pollard, 1988), which possibly created favorable stress conditions for silicic magmatic ascent and eruption at some locations (Bursik et al., 2003). Conversely, silicic magmatism can impose changes on the surrounding tectonic stress field, affecting fault behavior. Evolving magma pressures and crustal temperatures associated with the growth and development of large-scale silicic magma bodies will locally control crustal rheology and the resolved stress state in the crust and on surrounding faults (Simakin and Ghassemi, 2010; Gregg et al., 2012; Green et al., 2015; Ruz Ginouves et al., 2021). The subsequent triggering of fault slip on overlying caldera faults can in turn drive pressure changes in the underlying silicic magma system, with the potential to trigger voluminous silicic eruptions (Gudmundsson, 1998; Gregg et al., 2012). Large-scale caldera-forming eruptions also drastically alter the structure of the upper crust, driving large fault displacements triggered by the evacuation of large volumes of magma (Spray, 1997; Acocella, 2007; Gravley et al., 2007; Kennedy et al., 2018).

The complex tectonic-magmatic interactions that occur during large silicic eruptions demonstrate the critical need to understand the tectonic systems surrounding caldera volcanoes, as well as the history of the regional volcanic, tectonic, and volcano-tectonic events and the timescales over which they operate. Furthermore, given the tendency for frequent earthquake activity in regions undergoing extensional strain (e.g., Natron Basin, Tanzania; Lake Malawi Rift, Malawi; Long Valley, U.S.A.: Hill et al., 1995; Biggs et al., 2010; Kolawole et al., 2018; Reiss et al., 2021), it can be challenging to assess whether periods of apparent unrest are associated with magma and/or hydrothermal fluid migration, accommodation of tectonic extension, or a combination of these (Acocella et al., 2015; Illsley-Kemp et al., 2021; Reiss et al., 2021). This is particularly true of silicic caldera systems where the frequency of large earthquakes (M > 5) within close proximity to these volcanoes (i.e., <30 km) significantly exceeds that of volcanic eruptions (e.g., Long Valley Caldera, United States, and Corbetti volcano of the Main Ethiopian Rift: Hill et al., 2003; Martin-Jones et al., 2017; Wilks et al., 2017).

In this contribution we investigate volcano-tectonic interactions at silicic calderas by presenting an overview of deformation at Taupō, a large silicic caldera volcano in the central North Island of New Zealand (Figure 1). Comprehensive reviews of the nature of tectonic-magmatic interactions and deformation associated with caldera unrest are presented in a number of previous works (e.g., Newhall and Dzurisin, 1988; Gudmundsson, 2020; Acocella, 2021). Here, we review the patterns and potential sources of crustal-scale deformation at Taupō volcano using a wide range of observational tools (e.g., geological observations, historical reports, lake-leveling data, seismicity) to reveal how deformation can manifest on a range of timescales (annual to 10⁴ years) at a rifting silicic caldera, with much of this deformation driven by, or interacting with, the underlying magmatic system. Although there is no historically documented eruption from Taupō, this volcano has a well-constrained eruption history in its geological record (e.g., Barker et al., 2021) and numerous examples of recent seismic and surface deformation unrest events (e.g., Grindley and Hull, 1986; Illsley-Kemp et al., 2021). These data and observations allow us to examine the spectrum of volcanic-tectonic interactions at this silicic caldera to determine the characteristic deformation patterns and potential interactions between magmatic-tectonic processes during eruptive and non-eruptive unrest events.

FIGURE 1. Annotated DEM map of part of New Zealand’s North Island (inset bottom right) showing the generalized structure of the Taupō Volcanic Zone and other features discussed in the main text. Modeled geodetic extension rates (black arrows) for the onshore TVZ (in mm/yr) are from Wallace et al., 2004). Active faults identified and presented in the active fault database of Langridge et al. (2016) are represented by thin black lines. The position and structure of silicic calderas are generalized and based on Villamor et al. (2017). A more precise representation of the modern Taupō caldera is presented in Figure 2. TFB = Taupō Fault Belt.

2 Present Day State of Taupō Volcano

2.1 Geological Background

Taupō is one of two active rhyolitic caldera volcanoes in the Taupō Volcanic Zone (TVZ). The TVZ represents a complex region of volcanism and rifting that is related to oblique plate convergence and subduction of the Pacific Plate under continental New Zealand (Cole, 1979; Cole, 1990; Wilson et al., 1995) (Figure 1). To the north of New Zealand, in the Kermadec arc, volcanism is expressed through a typical arc front and back-arc system (Wysoczanski et al., 2010). However, the TVZ does not have a clear arc front to back-arc configuration and has instead been described as a “rifted arc” (Wilson et al., 1995; Wilson et al., 2009). The Hikurangi subduction zone, which drives the system, has an increase in trench convergence rates and steepening of the subducting slab to the north (Seebeck et al., 2014). This leads to clockwise rotation of the forearc block between the TVZ and the trench and rifting rates in the TVZ that change from near zero at the southern termination near Whanganui to up to ∼13–15 mm/yr near the Bay of Plenty coastline and offshore Whakatane Graben (Villamor et al., 2001; Wallace et al., 2004; Lamarche et al., 2006; Berryman et al., 2008; Gómez-Vasconcelos et al., 2017) (Figure 1). The TVZ is segmented by volcano type along its length, with andesite-dacite stratovolcanoes in the southern (e.g., Ruapehu and Tongariro) and northern (e.g., Whakaari) segments, and a central segment dominated by silicic caldera volcanoes (Wilson et al., 1995). The central TVZ, which includes Taupō volcano, is particularly notable in being the most active and productive region of silicic volcanism on Earth during the Quaternary (Wilson et al., 2009). The combination of rifting and enhanced fluid flux from the subducting slab plays a major role in driving the extreme heat flow and volcanism in this region (Barker et al., 2021, and references therein), and regions of active rhyolite volcanism are associated with higher heat fluxes and mantle melt production (Barker et al., 2020).

Along-axis variations in magmatism have been used to interpret the spatial and temporal variability of rifting across the broader TVZ (Rowland et al., 2010). Rifting in the central TVZ in particular is accommodated through a wide range of mechanisms including pure tectonic faulting, dike intrusion, and co-volcanic and/or co-intrusive fault slip, with regional scale uplift or subsidence due to the interplay between magma and rifting (Beanland et al., 1990; Rowland and Sibson, 2001; Gravley et al., 2007; Berryman et al., 2008; Rowland et al., 2010; Villamor et al., 2011; Villamor et al., 2017; Gómez-Vasconcelos, et al., 2020b). Volcano-tectonic interactions clearly play an important role in the style of rifting; however, many uncertainties remain over how these interactions vary over short (years), medium (10²–10³ years) and long timescales (10⁴–10⁵ years), particularly in response to the evolution of large silicic magmatic systems or changes in regional tectonic patterns. In addition, there is growing evidence that tectonics can play a major role in controlling the timing and behavior of eruptions and unrest for caldera volcanoes in this setting (Allan et al., 2012; Barker et al., 2016; Allan et al., 2017; Illsley-Kemp et al., 2019; Illsley-Kemp et al., 2021). To build a more complete picture of the dynamic magmatic-tectonic interactions occurring at rifted silicic calderas, we here integrate existing geophysical and geological data from Taupō volcano to provide a detailed picture of the volcano-tectonic interactions occurring at this system on annual to 10⁴ year timescales.

2.2 Surface Volcanics and Shallow Structures of Taupō Volcano

The general structure of Taupō volcano and surrounding regions reflects a history of volcanic-tectonic processes and interactions in this rifted arc setting (Figure 2). Rhyolite lavas and minor ignimbrites dated between ∼100 and >300 ka from Taupō occur over a wide region around the lake (Leonard et al., 2010; Barker et al., 2021). However, most of the young activity has been focused in the area that is defined by and concealed beneath Lake Taupō (Barker et al., 2021). Although Taupō volcano is situated at the southern end of the ∼340 ka Whakamaru caldera (Wilson et al., 1995), the modern caldera structure at Lake Taupō can be largely attributed to the 25,580 ± 258 years BP (rounded to 25.5 ka hereafter) Oruanui eruption, which discharged ∼530 km³ DRE (dense-rock equivalent) of magma in a complex 10-phase phreatoplinian eruption yielding widespread fall deposits and ignimbrite, together with thick (2–3 km) intracaldera fill (Wilson, 2001; Van Eaton et al., 2013). Since the Oruanui eruption, activity has been characterized by three small volume (<0.1 km³) dacites (units Ψ, Ω and A, from ∼20.5 to 17 ka) from vents at the northern end of the lake, and then 25 rhyolite eruptions between ∼12 and 1.8 ka, that range in volume over several orders of magnitude (units B through Z: Wilson, 1993; Sutton et al., 1995; Sutton et al., 2000; Barker et al., 2015; Barker et al., 2019; Barker et al., 2021). No relationship is observed between eruptive volume and repose period for these young eruptions. The most recent explosive eruption at 232 ± 10 CE (Taupō eruption, Unit Y: ∼35 km³ DRE; Wilson, 1993; Hogg et al., 2012; Hogg et al., 2019) represents the largest post-Oruanui event and resulted in further caldera collapse. More than 60% of the observed and inferred rhyolitic vent sites form a NE-trending lineament along the eastern side of the lake (Wilson, 1993; Barker et al., 2015; Barker et al., 2019) (Figure 2). This lineament is parallel to the general NE-trend of the TVZ rift and occurs at the eastern edge of the Oruanui caldera collapse structure identified from gravity data (Davy and Caldwell, 1998) (Figure 2). Almost all post-25.5 ka vent sites are contained within the Oruanui caldera, although small lava domes occur on the eastern margin of the lake and to the north of the lake in Acacia Bay (Figure 2). Heat flow data (Whiteford, 1996) and observed gas venting on the lake floor of Taupō (De Ronde et al., 2002) suggest an active hydrothermal system; however, very little is known about this system in part due to restricted access to the lake bed.

FIGURE 2. DEM hillshade and bathymetry map of the Lake Taupō region illustrating important volcano-tectonic features discussed in the main text. Orange shaded region represents the inferred footprint of the modern magma system from Illsley-Kemp et al. (2021). Caldera boundaries and post-Oruanui vent locations (orange-filled circles) are from Barker et al. (2021), and active faults are from the active fault database of Langridge et al. (2016). An intense negative gravity anomaly defines the Oruanui structural caldera (Davy and Caldwell, 1998), which is differentiated by the Oruanui collapse collar that extends to the eastern and western edges of northern Lake Taupō (Lipman, 1997). Black-filled circles represent the locations of rhyolitic domes of the NE dome magmatic system discussed in the text (e.g., RTD = Rubbish Tip Dome) and also the location of K-Trig basalts. Major faults of the Taupō Fault Belt are also annotated.

A distinctive TVZ rift fault pattern is expressed along the northern shore and north of Lake Taupō within the Taupō Fault Belt (Figures 1, 2). This region of faulting forms part of the Whakamaru subdomain of Rowland and Sibson (2001), which comprises a dominantly NE-trending structural fabric represented by normal fault traces, systematic fractures, and linearly aligned basaltic vents (Barker et al., 2021, and references therein). The major NE-trending faults controlling the structural relief are, from southeast to northwest, the Kaiapo, Whakaipo, Ngangiho, and Whangamata faults (Figure 2). Each of these faults has been associated with confirmed or inferred ground cracking during seismic unrest events documented during the 20th century (Hull and Grindley, 1984; Otway, 1986; Otway et al., 1984; Potter et al., 2015; see also Section 5). A subordinate N-trending set of faults is observed in the northwest (Figure 2). These faults have trends that parallel faults farther north in the Hauraki Rift (Figure 1), occur near the inferred outer edge of the Whakamaru caldera, and also align with a N-trending elongate cluster of rhyolite domes (Rowland and Sibson, 2001). However, the N-trending faults have also been defined as inactive based on the current definition of fault activity for TVZ faults in the New Zealand active fault database (Langridge et al., 2016) and, as they do not displace Oruanui deposits, likely have not ruptured in the last ca. 25,500 years (Villamor and Berryman, 2001).

The expression of Lake Taupō south of the Oruanui structural caldera resembles a typical linear rift lake (e.g., Lake Tanganyika; Wright et al., 2020; Shaban et al., 2021). Subsidence in southern Lake Taupō occurs largely on the downthrown side of the Waihi fault, with the footwall of this structure forming the NE-trending Karangahape Cliffs on the western edge of the lake (Wilson, 2001) (Figure 2). Similarly, a series of NE-trending faults control the eastern expression of the Lake Taupō shoreline. Although the structural expression of Lake Taupō in many places has the trademarks of a NE-trending rift system, much of this structural relief, particularly since the Oruanui eruption, is inferred to have developed in response to caldera collapse (e.g., Davy, 1993; Wilson, 2001; Barker et al., 2021).

Previous studies suggest that the tectonic structures surrounding Lake Taupō can play a critical role in events leading to silicic eruptions (Allan et al., 2012; Barker et al., 2016). For example, an active magmatic system (NE dome magmatic system; Figure 2) defined by elevated heat flux within active geothermal systems (Barker et al., 2021; Bibby et al., 1995; Rosenberg et al., 2020) likely occurs to the northeast of the current magma mush system below Lake Taupō (Barker et al., 2015), but has been relatively quiet since the Oruanui eruption with only one eruption at ∼16.5 ka (Puketarata; Kósik et al., 2019; Figure 2). However, lateral movement of shallow magma from this system towards Taupō parallel to the NE-trending structural fabric is inferred to have occurred during the Oruanui eruption (Allan et al., 2012), suggesting that dike-driven rifting events can be associated with voluminous silicic eruptions at Taupō, while revealing a local interplay between different magmatic and tectonics systems at >10 km length scales (cf. Biggs et al., 2016).

2.3 Taupō Crustal Magmatic System

Geochemical studies of Taupō volcano’s eruption products provide critical constraints on the subsurface structure in the <10 km-deep shallow portion of the crust, distribution of partial melt, and recent state of the shallow magma system (e.g., Barker et al., 2015; Barker et al., 2016) (Figure 2). Since the Oruanui eruption, a large silicic magmatic system has rebuilt at Taupō and discharged three closely related rhyolite magma subgroups in multiple eruptions between ∼12 and 1.8 ka (Barker et al., 2015). Mineral barometry indicates crystallization at a range of pressures for rhyolites of the first subgroup, erupted from 12 to 10 ka, but dominantly at pressures of ∼100–200 MPa (∼4–8 km depth) for the second and third subgroups. These depth estimates, along with the superimposed vent and collapse areas highlight that the young system overlaps in its subsurface extent with the former Oruanui magmatic system (Bégué et al., 2014; Barker et al., 2015; Myers et al., 2018; Myers et al., 2019). The range of crystallization depths for the post-Oruanui Taupō magmatic system, combined with the areal distribution of vent locations, support estimates of a maximum volume of mush within the present system of 200–1,000 km³ (Barker et al., 2015). The spatial extent of the modern magma system has also been estimated through seismic analysis of the 2019 unrest episode (described in detail in Illsley-Kemp et al., 2021) (Figure 2), where a distinct gap in seismicity occurred over an area of >80 km². The seismic data lead to an inference that the current shallow silicic magmatic system likely contains >20–30% melt and has a minimum volume of 250 km³ (Illsley-Kemp et al., 2021). The shallow silicic system is fed by a deeper reservoir and mafic magma recharge from depth, with juvenile mafic enclaves and basaltic recharge evident in larger past eruption products (Blake et al., 1992; Barker et al., 2016, 2020).

2.4 Deep Crust and Uppermost Mantle Below Taupō

The magmas feeding the Taupō magmatic system are sourced from the Hikurangi mantle wedge. This region is characterized by high seismic attenuation and low seismic velocities, consistent with the presence of partial melt and fluids from the subducting Hikurangi Plateau (Reyners et al., 2006; Eberhart-Phillips et al., 2008; Eberhart-Phillips et al., 2013; Eberhart-Phillips et al., 2020; Eberhart-Phillips and Fry, 2018). The compositions of primitive magmas reaching the surface highlight that both rift-induced decompression melting and flux melting play a role in TVZ magma genesis (Barker et al., 2020), with higher degrees of partial melting and higher magma flux rates into the crust in areas where silicic caldera volcanism dominates.

Geophysical monitoring stations (e.g., seismometers) are restricted to onshore Lake Taupō, and thus the resolution and extent of geophysical imaging immediately below Lake Taupō is limited. Our understanding of the lithospheric structure and distribution of crustal and upper mantle melts is instead largely informed by geophysical surveys immediately north of the lake (e.g., Bannister et al., 2007; Stern and Benson, 2011), supplemented by seismicity data and geochemical and petrological observations of erupted products (e.g., Barker et al., 2015; Illsley-Kemp et al., 2021) (Figure 3). Analysis of active source data has revealed a zone of increased P-wave speeds (>6.8 km/s) below the TVZ at 15–20 km depth and V_p exceeding 7.3 km/s below 20 km, which was interpreted by Stratford and Stern (2006) to reflect the presence of anomalous mantle. However, an alternative model was developed by Harrison and White (2006) for the same area, based on active and passive source seismic data, where the upper mantle was interpreted to start at 30–35 km depth and the high P-wave speeds between 15 and 30 km were thought to represent the presence of dense magmatic underplate material.

FIGURE 3. The geophysical structure of the crust beneath the Taupō Fault Belt, adapted from Stern and Benson (2011) (position of cross-section line shown in Figure 1). The rift pillow has anomalously low seismic velocities and is interpreted to represent a region of partially molten mafic intrusions. The mantle reflector is thought to be a region of highly molten primitive basalt. Topography is 10x exaggerated.

To resolve this discrepancy, Stern and Benson (2011) adopted wide-angle seismic reflection/refraction methods to image in greater detail the crust and upper mantle structure along an NW-SE transect ∼5 km north of Lake Taupō (Figures 1, 3). From analysis of reflected events they interpreted the lower crustal and upper mantle structure in the vicinity of Taupō volcano, and the distribution of high temperature zones, fluids and partial melt from their seismic velocity estimates. By considering their results in the context of existing magnetotelluric data (Heise et al., 2007) and other independent seismic studies (Bannister et al., 2004; Behr et al., 2011), their key findings included: 1) a high V_p region of variably cooled mafic material between 15–25 km depth, which occurs immediately below a low resistivity zone of pervasive partial melt at 10–15 km depth (Heise et al., 2010); 2) a boundary to the upper mantle with an unusually low seismic velocity at a depth of 25 ± 2 km beneath the rift axis; and 3) a highly reflective zone at ∼30 km depth interpreted to represent the upper contact of a lens of pooled basaltic melt in the uppermost mantle (Figure 3).

3 Thirty-Five Years of Taupō Seismicity

Earthquake patterns in the Lake Taupō region reflect an ongoing interplay between magmatic, tectonic, and geothermal systems (Grindley and Hull, 1986; Potter et al., 2015; Illsley-Kemp et al., 2021). Since the mid-1980s, New Zealand’s permanent seismic monitoring network has allowed for continuous earthquake monitoring around Lake Taupō, covering seven unrest episodes that exhibited enhanced seismicity (Figure 4). Large magnitude earthquakes showing classic foreshock-aftershock sequences delineating fault-planes are rare, and seismicity more often occurs in clustered swarms dominated by relatively low magnitude events, indicative of fluid-driven fracturing (e.g., Lindenfeld et al., 2012; Jacobs et al., 2013; Bannister et al., 2016; Reiss et al., 2021) associated with the movement of either magma and/or hydrothermal fluids. This interpretation is further supported by full moment tensors calculated in a series of swarms below Lake Taupō in 2019, where the majority of events exhibit non-double couple solutions (Illsley-Kemp et al., 2021).

FIGURE 4. Locations and timing of seismic swarms in the Lake Taupō region from 1985, when the local seismic network was improved, to June 2021. (Panels A to G) show the locations of earthquakes in each swarm, and (Panel H) shows all other earthquakes not associated with these swarms (i.e., background seismicity). All earthquake information is taken from GeoNet (Petersen et al., 2011), the magnitude range is 2–5 and the depth range is 0–25 km. (Panel I) shows the cumulative number of earthquakes and seismic moment release through time, with the duration of each swarm highlighted in blue.

Generally, seismicity below the lake occurs outside the inferred boundary of the Oruanui structural caldera (Figure 4G) and exhibits recurring patterns, such as a roughly WNW-ESE trending cluster near the southern margin of the caldera, as well as seismic clustering around the NE-SW trending Horomatangi Reefs (Figures 4D,E,G) and in the vicinity of the youngest vent sites for the volcano (Barker et al., 2021). Additional swarms have occurred outside the lake area to both the north and south. Seismicity north of the lake from December 2000 to June 2001 occurred within the Taupō Fault Belt (Figure 4C), during which time uplift was detected immediately east of the Kaiapo Fault (within ∼5 km; see also Section 5 and Figure 5). Swarms south of the lake are characterized by a high number of earthquakes typically exhibiting low magnitudes (Figures 4E,F). The largest magnitude earthquakes experienced around Lake Taupō since 1985 include M5.2 and 4.8 events in 2019 (Illsley-Kemp et al., 2021), and M4.2 and 4.6 earthquakes in 2002 that were not associated with swarm activity.

FIGURE 5. Relative elevation changes recorded on lake leveling stations around Lake Taupō [elevation changes are relative to station WA in (Panel C)]. Each colored line represents relative changes in elevation through time for an individual station, with the specific station indicated by the symbology on the right side of the graphs in (Panels A and B). Black dashed line in each plot signifies the initiation of the 1983 Kaiapo seismic swarms. The stations are divided into (A) “Fault Belt” sites, which represent stations north of the lake that are positioned within the Taupō Fault Belt, and (B) “Near Vent” sites, which represent stations in close proximity to the eruptive vent locations for the 232 CE Taupō eruption and the Horomatangi Reefs (Figures 1, 2). Station RA is colored both red and blue as it occurs within the Taupō Fault Belt and near to the Horomatangi Reefs. The locations of each station relative to the shoreline of Lake Taupō are shown in (Panel C). All lake-leveling data used in this study can be accessed from DOI: 10.5281/zenodo.5777476.

In rare instances, such as 2000–2001, seismicity exhibits a mainshock-aftershock sequence and appears to have a purely tectonic origin (Utsu and Ogata, 1995), while other seismic swarms do not show any characteristic tectonic signature (Figure 4; see also Jacobs et al., 2013). Activity during month-to year-long seismic unrest periods manifests as a number of individual day-to week-long seismic swarms that define distinct geographical clusters (e.g., Figure 4G). In detail, seismic swarms can exhibit rift-parallel sub-clusters that are ∼2 km wide and up to ∼6 km long (Illsley-Kemp et al., 2021). These linear arrays typically trend NE-SW and occasionally NNW-SSE; the NE-SW trending swarms can occasionally be linked to faults observed at the surface, supporting the active role of rift-aligned faults in accommodating regional tectonic and/or local magmatic stresses. The volcanic unrest event in 2019 also caused seismicity and ground deformation centered around the Taupō collapse caldera, with most seismicity occurring near the southern edge of the interpreted magma reservoir and exhibiting a rough ENE-WSW trend (Illsley-Kemp et al., 2021). In all, these data support the role of both tectonic and magmatic processes in driving the ongoing deformation in the region.

4 Surface Deformation Patterns From Taupō Lake Leveling

The longest record of detailed ground deformation data from Taupō is from lake leveling surveys, which began in April 1979 (e.g., Otway, 1980; Otway, 1986; Otway et al., 2002). The data are gathered using a portable gauge at 22 localities around the lake shore (Figure 5) and corrections account for lake seiching and variable water levels. Although previously published as incremental updates, here we present the dataset in its entirety with all offsets given relative to site Waihaha (WA) in the Western Bay area of the lake (Figure 5), which occurs to the west of the Taupō Fault Belt and outside the boundaries of modern rifting delineated by Villamor et al. (2017). This allows us to compare vertical deformation between regions within and outside of the rift. Note that previous descriptions of the dataset (e.g., Peltier et al., 2009) have used Rainbow Point (RP) as a reference site; however, it is likely that this locality is itself affected by regional deformation (Otway et al., 2002). For simplicity we select lake leveling localities in two groups; those that are north of the lake in the Taupō Fault Belt (“Fault Belt” sites in Figure 5A, red boxes), and those that are close to the Horomatangi Reefs (“Near-Vent” sites in Figure 5B, blue boxes).

Since the lake leveling surveys began, the first major deformation was observed in March 1983, when an uplift rate of >5 mm per month was observed at Kinloch (KH), which prompted the installation of new lake leveling localities. After an earthquake swarm in June 1983 (Webb et al., 1986), which was followed 1 week later by ground rupturing along the Kaiapo Fault (Hull and Grindley, 1984), Kinloch (KH) and the neighboring sites (WO, TA) began subsiding. Meanwhile Rangatira (RA) rose steadily after October 1983, peaking at 60 mm uplift by January 1984 (Figure 5B).

After the events of June 1983, all sites NW of the Kaiapo Fault (downthrown side) have subsided on the order of millimeters per year at near-constant rates, although localities NW of the Whakaipo Fault (e.g., WO, TA, KH) have subsided at faster rates (up to 7 mm/yr, compared to ∼2 mm/yr at KO and WI), suggesting that both the Kaiapo and Whakaipo Faults were possibly accommodating some deformation during this time. In June 2000, sites RA and AB east of the Kaiapo Fault began rising, followed in 2001 by a ∼1 year period of uplift across all the Fault Belt sites, and coincident with a major earthquake swarm just north of the lake (Figures 4C, 5A). This uplift continued for all sites SE of the Whakaipo Fault until late-2008, when it then reversed. This region has remained stable since then. The uplift period between 2001 and 2008 at some Fault Belt sites is approximately coincident with uplift near the Horomatangi Reefs (Figures 5A,B) suggesting a possible volcanic source, with the termination of uplift in late 2008 also coinciding with heightened seismic activity within the caldera (Figure 4D).

Vertical elevation changes in Fault Belt sites over the last 40 years show a clear bimodal pattern either side of the Kaiapo Fault (fault location shown in Figure 2), with uplift to the east and subsidence to the west of the structure. However, the overall data show no evidence for a discrete dislocation along the Kaiapo Fault plane, with the relative displacement showing a much longer wavelength than those characteristic of discrete fault slip events (e.g., Jara-Muñoz et al., 2019). Although the Kaiapo Fault represents the largest surficial fault scarp in the Taupō Fault Belt and has documented surface breaks in historical times, the majority of vertical displacements in the last 40 years in North Taupō cannot be attributed solely to slip along this fault based on the presented lake levelling data.

5 Historic Unrest and Surface Rupturing Events at Taupō

The Lake Taupō region exhibits a history of persistent earthquakes and related surface deformation (Grindley and Hull, 1986), suggesting that the area is prone to earthquake hazards while also highlighting the potential role of tectonic triggering for past eruptive events (see also Allan et al., 2012; Barker et al., 2016). Sixteen periods of unrest have been identified by Potter et al. (2015) as occurring between 1872 and 2011. Here, we discuss only the most intense periods of recorded faulting and associated ground shaking in regions up to ∼15 km from the edge of the Taupō caldera. Three of these events resulted in observed surface rupturing, in 1895, 1922 and 1983, whereas no surface breaks were recorded in association with seismicity and ground shaking during 1964–1965 and 2019.

Seismicity in 1895 lasted for 6 weeks, producing ground cracks, landslides and subsidence in the broader Lake Taupō region. The largest earthquake on August 17 was also felt in Napier, Hamilton, and Gisborne, with reported damage indicating shaking intensity of MM8 (Modified Mercalli scale) and an earthquake magnitude likely exceeding M6 (Eiby, 1968). In 1922, seismic activity was recorded from May to October, and was accompanied by widespread surface deformation along a number of known faults in the region (Morgan, 1923; Grange, 1932) (Figure 6). The largest earthquake in the sequence on June 10 was associated with surface displacements on the Kaiapo Fault, and surface ruptures as long as 3 km were observed on the Kaiapo, Whakaipo, and Ngangiho Faults, with ground cracking observed over a ∼140 km² area (Figure 6). Recorded vertical displacements typically ranged between 0.1 and 0.9 m, with displacements as high as 3 m reported at some localities. However, the extent to which these vertical displacements record pure fault-slip activity vs. other near-surface processes (e.g., secondary landsliding, sediment compaction) is unknown, particularly for larger documented displacements.

FIGURE 6. Historical surface deformation associated with the 1922 (Panel A) and 1983 (Panels B and C) seismic swarms. Distributions of fault ruptures, mapped surface cracks, locations of reported (but not mapped) surface cracks, and regions of reported subsidence are summarized and documented in Morgan (1923); Grange (1932); Hull and Grindley (1984), and Villamor et al. (2001). Images in (Panels B and C) are from Otway et al. (1984) and show surface ruptures that formed along the Kaiapo Fault with both cm-scale opening and vertical slip components during the 1983 seismic event.

Activity in 1983 occurred over an extended period, with two main seismic swarms recorded in February and June-July. The “Kaiapo swarm” in June-July was associated with the two largest earthquakes of M ∼4.3, and 64 events of M > 3. On July 1, a ∼1 km deep, M 4.3 earthquake accompanied rupture of the Kaiapo Fault at the same location recorded in 1922, producing ∼1.2 km of en-echelon surface fracturing with 30 mm maximum extension and 50 mm maximum throw (Otway et al., 1984; Grindley and Hull, 1986) (Figures 6B,C).

5.1 Subsidence and Uplift Patterns During Historic Unrest Periods

Periods of seismic unrest at Taupō are occasionally associated with detectable surface uplift and/or subsidence, which likely indicate ongoing magmatic, hydrothermal, and/or tectonic activity in the region (Peltier et al., 2009). Ground subsidence was observed in the Whakaipo Bay region for both the 1922 and 1983 seismic swarms (e.g., Grindley and Hull, 1986) (Figure 6A). In 1922, an apparent lake level drop was reported prior to the main seismicity event in June, likely indicating decimeter-scale uplift. However, by the end of June of 1922, the Whakaipo Bay region had subsided 0.6–0.8 m, and the region had subsided 2.3 m on average by the end of 1922. A few years later, Grange (1932) measured a total subsidence of 3.7 m in the Whakaipo Bay region since the initial major period of seismicity in June 1922.

A similar uplift and subsidence pattern is reported for the 1983 seismic event, although the scale of this deformation is on the order of centimeters rather than meters. Lake levelling data reveal up to 43 mm of uplift in North Taupō from 25 October 1982 to 4 June1,983, with maximum uplift centered in the Kinloch region ∼10 km west of the Kaiapo Fault (Otway, 1986). As much as 14 mm of subsidence was also observed at this time at the southern end of the lake. These surface deformation patterns changed during and after the 1983 June-July seismic swarm (Grindley and Hull, 1986; Otway, 1986). Lake leveling data from June 4 to 19 July 1983, revealed up to 55 mm of subsidence on the hanging wall side of the Kaiapo Fault and 10 mm of uplift on the footwall, consistent with the measured throw values observed on surface ruptures (50 mm maximum throw: Otway et al., 1984).

6 Long-Term Deformation at Taupō Volcano

6.1 Time-Average Fault Slip Rates Along the Taupō Fault Belt

Large eruptions associated with ignimbrite formation, such as the 25.5 ka Oruanui eruption, tend to overwhelm the landscape and overprint and erase the earlier faulted topography. As such, minimum time-averaged fault slip rates since ∼25.5 ka can be estimated by measuring the vertical offsets of fault scarps dissecting the Oruanui ignimbrite surface (Figure 7). Across-strike real-time kinematic GPS transects measuring surface fault displacements along the Taupō Fault Belt presented in this study (Figure 7, see also Supplementary Material) support time-averaged vertical slip rates of 0.6 ± 0.07, 1.1 ± 0.13, 1.2 ± 0.11, 0.9 ± 0.09, and 0.9 ± 0.09 mm/yr for the Aratiatia, Kaiapo, Whakaiapo, Ngangiho, and Whangamata faults, respectively, accounting for a cumulative vertical slip rate of 4.6 ± 0.23 mm/yr (Figure 7B). Assuming fault dips of 60°, these cumulative fault displacements across the Taupō Fault Belt correspond to 2.7 ± 1.2 mm/yr extension in the shallow crust along a ∼15 km-wide zone of faulting. Adopting the approach of Villamor and Berryman (2001) (see also Supplementary Material), these surficial measured extension rates could account for as much as 5.7 ± 2.6 mm/yr extension within the seismogenic crust, depending on how fault dips and the total fault slip change with depth. These results suggest that a significant portion of the ∼10 mm/yr of extension implied from GPS data (Wallace et al., 2004) is accommodated immediately north of Taupō along the Taupō Fault Belt, with likely a subordinate proportion of this extension also accommodated in the present day in the Taupo-Reporoa Basin (Figure 1). However, it is unclear what proportion of this near-surface extension observed along faults might be accommodated through dike intrusion at depths >1 km (i.e., Acocella and Trippanera, 2016), as dike-induced faulting represents a common mode of extension accommodation for magmatic rift systems elsewhere (e.g., Iceland, Afar Depression (Ethiopia), Kenya Rift (Kenya); Rubin and Pollard, 1988; Rowland et al., 2007; Muirhead et al., 2016).

FIGURE 7. Long-term deformation trends around Lake Taupō and Taupō Fault Belt. (A) DEM hillshade map showing the position of the fault throw transect, which involved measuring the heights of observable fault scarps (1 m minimum throw) displacing 25.5 ka Oruanui eruption deposits using real-time kinematic GPS. The locations of the active faults measured here, and shown on the map, are from the active fault database of Langridge et al. (2016). Note that, at the location of the transect, the primary faults of the Taupō Fault Belt (e.g., Whakaipo and Kaiapo faults) are composite faults with a number of splay structures. (B) Cumulative vertical fault slip rate measured on all fault scarps shown in the transect in Panel A displacing the Oruanui surface. These vertical slip rate values are time-averaged, based on the surficial slip accrued on each fault since emplacement of the Oruanui deposits. The cumulative vertical slip rate has an estimated error of 45% (see Supplementary Material). (C) Cumulative uplift and subsidence rates along the transect shown in Panel A, based on fault slip measurements on the Oruanui surface. All rates are relative to the footwall of the northwest-most fault. Results reveal a zone of greatest relative subsidence along the rift axis. (D) Map showing the locations of the recorded elevations of the post-Oruanui shoreline highstand from Manville and Wilson (2003). Also shown for reference are the primary faults of the Taupō Fault Belt. Compass wheel relates to the bearings shown in Panel E. (E) Elevations of the post-Oruanui shoreline highstand and the positions of faults are from Manville and Wilson (2003). The positions of these features relate to the bearing from the center of the lake (see also compass bearings shown in Panel D).

6.2 Surface Deformation From Paleoshorelines

The two major silicic eruptions, the 25.5 ka Oruanui (530 km³ DRE) and 232 CE Taupō (35 km³ DRE) eruptions, respectively, at Taupō volcano were followed by significant changes to paleohydrologic systems surrounding the lake (i.e., pre-Oruanui “Lake Huka” and modern-day Lake Taupō: Wilson et al., 1997; Manville et al., 1999; Manville and Wilson, 2003; Barker et al., 2021). Surface drainage networks were altered or destroyed as pyroclastic materials filled and dammed nearby valleys and channels (Smith, 1991), ultimately increasing the post-eruption spill elevation of the lake (Wilson et al., 1997; Manville and Wilson, 2004). These post-eruptive rises in lake level are indicated by variably developed terraces and associated beach deposits around the lake, revealing that the paleoshoreline of Lake Taupō was situated as high as ∼140 m and ∼34 m above the modern-day lake level in the immediate aftermath of the Oruanui and Taupō eruptions, respectively (Manville and Wilson, 2003; Manville and Wilson, 2004). However, within decades of the eruptions, barriers damming the Lake Taupō outlet (i.e., the Waikato River) were breached and then rapidly degraded and destroyed (Manville et al., 1999; Manville and Wilson, 2004), restoring the lake level closer to its pre-eruptive elevation. Therefore, due to the close temporal associations between high-stand paleoshoreline development and the major eruptions at Lake Taupō, paleoshorelines represent excellent time markers for examination of time-averaged surface deformation patterns since the Oruanui and Taupō eruptive events (Figures 7D,E). Here we focus specifically on the post-Oruanui paleoshoreline to examine how caldera magmatism may impact nearby rifting processes on 10⁴ year timescales.

GPS-based full crustal extension rates in the Lake Taupō region (Wallace et al., 2004) and time-averaged extension rates estimated from fault displacements on the Oruanui ignimbrite support between 10 and 5 mm/yr, respectively, of extension immediately north of Lake Taupō (Figures 7A–C). Measured time-average vertical surface displacements in the Taupō Fault Belt presented here support cumulative vertical slip rates of ∼4–5 mm/yr at the surface (Figure 7B), which over thousands of years should produce subsidence along the rift axis on the order of a few mm/yr (Rowland and Sibson, 2001; Villamor and Berryman, 2001) (Figure 7C). However, the relative displacements of the post-Oruanui paleoshoreline markers around Lake Taupō do not show this expected pattern of rift axis subsidence (Manville and Wilson, 2003). For example, ∼25,500 years of axial rifting and related subsidence along the Taupō Fault Belt (minimum expected subsidence of 1–2 mm/yr from rift extension: Figure 7) would produce ∼25–50 m of relative vertical displacement and rift subsidence. Instead, elevations of post-Oruanui paleoshorelines reveal no appreciable net subsidence along the rift axis since ∼25.5 ka (Figure 7E), suggesting that the subsidence signature expected from rift-related faulting and crustal thinning has been overprinted due to magma emplacement into the crust and related surface uplifts (Manville and Wilson, 2003). These interpretations are supported by decadal-scale lake leveling observations (Figure 5; Section 4), suggesting that regions of Lake Taupō transition between phases of uplift and subsidence, likely resulting from a combination of rift extension, magma cooling and contraction, and variations in magma flux into the crust. Of note, these long-term uplift and subsidence patterns discussed above represent only those immediately north (within 10 km) of the Taupō caldera; within the caldera itself, caldera collapse resulting from the expulsion of >580 km³ of crustal material at and since ∼25.5 ka represents the primary control on the long-term subsidence pattern.

7 Syn-Eruptive Deformation

The presented observations (Sections 3–6) strongly suggest that episodes of surface deformation, seismicity and faulting in and around Taupō volcano are frequently modulated, and initiated, by subcrustal magma and fluid migration (Grindley and Hull, 1986; Smith et al., 2007; Peltier et al., 2009; Illsley-Kemp et al., 2021). In addition, a number of lines of evidence support significant syn-eruptive crustal and surface deformation during silicic eruptions, with deformation extending well beyond the caldera boundary, which we discuss below.

Tephra-filled fissures observed in exposures outside of the caldera represent evidence for strong ground shaking and related surface deformation during the 232 CE Taupō eruption. These fissures are primarily observed within local pre-Holocene loess (overlain by the ∼10 ka Unit E [Opepe] and ∼3.4 ka Unit S [Waimihia] deposits) and contain sharp straight to curviplanar matching walls that are commonly lined with a few millimeters thick layer of fine ash (Figures 8B,C). Infilling these fissures is the Taupō ignimbrite (Figure 8B), still preserving its original primary stratigraphy (i.e., layers 1 and 2: Wilson, 1985; Hughes, 2005). Similar Oruanui tephra-filled fissures are also observed 15 km east of Mt. Ruapehu in the upper Waikato Stream (Gómez-Vasconcelos et al., 2016) (Figure 8C). The field relationships discussed above suggest that these fissures represent vigorous ground shaking for distances up to 50 km from Taupō caldera, and were opened during or immediately prior to the Taupō and Oruanui eruptions and then progressively filled during these events.

FIGURE 8. Syn-eruptive deformation features associated with 25.5 ka Oruanui and 232 CE Taupō eruptions. (A) DEM hillshade map showing the position of modern Lake Taupō (black polygon) and the deformation features shown in Panels B–D. (B) Dilational fracture filled by the Taupō ignimbrite. The excellent wall preservation in a relatively unconsolidated host suggests the fracture opened shortly before and/or during emplacement of the ignimbrite. (C) 25.5 ka Oruanui tephra filling older tephra and lahar deposits (image courtesy of M.G. Gómez-Vasconcelos). As in Panel B, the fissure exhibits excellent wall preservation suggesting the fracture opened shortly before and/or during emplacement of the tephra. (D) Fault displacements accrued during and/or immediately after the Taupō eruption from Leonard et al. (2004) (image in the top panel, interpretation in the bottom panel). The N-S orientation of the outcrop face exposure is ∼45° to the dominant rift trend, and faults here represent a mixture of rift-parallel to rift-oblique structures (Leonard et al., 2004), suggesting the development of some faults may not be related exclusively to rift extension processes. Here, displacements on the Taupō ignimbrite do not extend into the overlying post-232 CE highstand lake deposits, suggesting that faulting occurred during and/or immediately after the eruption. The presence of a horst and graben sequence results in no net subsidence either side of the ∼35 m-wide fault system.

A region of intense syn-eruptive deformation (Figure 8D) is documented approximately 15 km NNE of the Horomatangi Reefs (Figure 1), aligning along-strike with the inferred positions of the 232 CE Taupō eruptive vents (Leonard et al., 2004; Houghton et al., 2010) (Figures 8A,D). Here a series of normal faults are observed in unconsolidated pyroclastic deposits (including those from all but the last phase of the 232 CE eruption) with displacements ranging from a few cm to 1.8 m (Figure 8D). These faults, in addition to some locally boudinaged indurated lacustrine sediments, accommodate ∼7 m of horizontal extension across the ∼35 m-wide outcrop (Leonard et al., 2004), but also show no net vertical displacement across the width of the local fault system. Across this exposure, equivalent fault displacements are observed on tephra units deposited prior to and during the 232 CE Taupō eruption, whereas the entire faulted sequence is truncated by post-Taupō beach and near-shore lacustrine sediments (Figure 8D), necessitating that this faulting developed either during or within years of the eruption (Leonard et al., 2004). The apparent lack of net vertical displacement suggests that the extension was driven, or additionally accommodated, by shallow magma emplacement (e.g., Magee et al., 2017) rather than pure tectonic faulting.

Collectively these field observations demonstrate that large silicic eruptions may be accompanied by significant extensional strain accommodation immediately outside of the caldera. Whether the fault structures we observe acted to accommodate regional tectonic extension or local magmatic stresses (i.e., dike-induced faulting: Rubin and Pollard, 1988), or both is, however, still largely unresolved and requires further investigation.

8 Discussion

8.1 Impact of Taupō Volcano on Rift Evolution and Behavior

The observations presented from Taupō provide key insights into how large silicic caldera volcanoes impact the evolution of continental rift architecture on annual to 10⁴ yr timescales. Ongoing extension at continental rift settings is accommodated primarily through normal faulting in the shallow subsurface (upper few kilometers), producing ongoing fault-driven subsidence that is either asymmetric and primarily focused on a single fault bounding a half-graben basin, or symmetric and focused within the rift axis (Morley, 1989; Ebinger, 2005; Muirhead et al., 2016). The subsurface stratigraphic position of syn-rift volcaniclastics in the central TVZ necessitates rift subsidence rates that are likely of the order of millimeters per year (Wilson et al., 2010; Downs et al., 2014; Villamor et al., 2017). Along the Ngakuru Graben and northern parts of the Taupō Fault Belt vertical fault displacements require at least 40 m of relative vertical subsidence focused along the rift axis within the last 20,000 years (Villamor and Berryman, 2001), and similar subsidence patterns are observed in faulted Oruanui ignimbrite deposits 10–15 km north of the lake (Figure 7). GPS and InSAR data also record recent short-term subsidence along the TVZ, although a proportion of subsidence observed today is attributed to cooling and contraction of axial magmas (Hamling et al., 2015), or alternatively deeper-seated changes in mantle flow viscosity (Lamb et al., 2017). The broad pattern of long-term rift-axis subsidence along the TVZ (Chambefort et al., 2014; Wilson et al., 2010) is not observed in the Taupō Fault Belt immediately north of Lake Taupō, where post-Oruanui paleoshorelines reveal no appreciable subsidence since ∼25 ka, and even show localized long-term uplift in places (Manville and Wilson, 2003; Figure 7). Additionally, the distribution of faulting in the Taupō Fault Belt immediately north of Lake Taupō (Figure 2) is more akin to block-style faulting (i.e., a series of horsts and grabens: Manville and Wilson, 2003), rather than an inward facing fault population (i.e., faults consistently dipping towards the central rift axis) that would drive axial subsidence.

Localized uplift around Taupō volcano is also recorded on decadal timescales. Almost 40 years of lake levelling measurements reveal continuous uplift at Near Vent sites centered around the Horomatangi Reefs in the east of the lake (Figure 5B). Net subsidence is recorded in regions northwest of the Kaiapo Fault, primarily resulting from deformation initiating sometime before the 1983 seismic unrest event (Hull and Grindley, 1984), which may have been caused by “dewatering” of an intrusion prior to 1979 or even 1963 (Smith et al., 2007). However, southeast of the Whakaipo Fault this subsidence ceased in 2001, reverting to the general uplift pattern seen in the northeast of the Lake Taupō shoreline. Overall, lake leveling data favor net uplift below and immediately north of the lake in the last 40 years, rather than the expected tectonic-driven subsidence due to rifting.

The described long-term deformation patterns are best explained by magmatic crustal additions (i.e., shallow sills and magma recharge into the mush system) within and around Taupō volcano’s magmatic system, which also help to maintain crustal thickness and limit rift subsidence (Manville and Wilson, 2003; Bastow and Keir, 2011). Thus, these deformation trends, which are not consistent with the expected longer-term (∼10⁴ yr) behavior of a continental rift (i.e., rift subsidence up to a few millimeters per year; e.g., Magadi Basin (Kenya) and Taupō-Reporoa Basin (TVZ); Wilson et al., 2010; Downs et al., 2014; Muirhead et al., 2016; Figure 7), are best explained by high crustal magma influx suppressing axial subsidence of the rift basin. Geochemical analysis of primitive basaltic magmas feeding the TVZ indicates that the degree of partial melting and mantle melt supply is 4–10 times higher beneath Taupō caldera when compared to the intercaldera regions farther to the north (e.g., Ngakuru Graben) (Barker et al., 2020) (Figure 1). Elevated magma flux therefore ultimately determines the location of the large caldera volcano, which in turn plays the primary role in controlling the structural expression of rifting within ∼10 km of the caldera (Figure 7). Between major silicic eruptions, magmatic intrusions produce localized uplifts that inhibit axial subsidence, while during eruptive events, meter-scale crustal extension can be accommodated north of the caldera along multiple faults with no appreciable subsidence observed across the fault system. This is illustrated by the syn-eruptive normal faults shown in Figure 8D, which form a horst and graben system with a cumulative vertical displacement greater than 5 m on the Rotongaio tephra, yet no vertical elevation change either side of the 35 m-wide fault system. Although localized subsidence is observed in between some faulted blocks, there is no net subsidence across the exposed fault system (see Figure 8D for illustration).

The observations discussed above illustrate that the impact of Taupō volcano on rift processes varies laterally with distance from the modern silicic magma system. What is not considered is how the areal extent of the silicic system may vary temporally, and its potential impact on crustal rheology and rift processes. Geophysical and geochemical observations constrain the maximum modern-day areal extent of the crustal magma system to a ∼250 km² region beneath Lake Taupō (Barker et al., 2021; Illsley-Kemp et al., 2021) (Figure 2). Furthermore, seismicity patterns suggest that heat exchange at the lateral edges of the magma mush affects the crustal rheology, with a transition from ductile to brittle behavior with distance from silicic mush (Illsley-Kemp et al., 2021). The accumulation of a large melt reservoir prior to eruptions at Taupō (Allan et al., 2017) could widen the areal extent of the current silicic system and will certainly increase the magnitude of heat exchange, resulting in transient rheological changes along the Taupō Fault Belt, which tentatively could affect the overall style of rifting immediately north (within 10 km) of the Taupō caldera boundary.

8.2 Volcanic-Tectonic Interactions at Taupō Volcano During Periods of Unrest, Eruption, and Non-Eruption

Rowland et al. (2010) proposed that the physical state of the mid-crustal hot zone in the central TVZ exerts a profound influence on rifting mechanisms including tectonic faulting, dike intrusion, magma body inflation and caldera collapse. These mechanisms are considered to operate on cycles that are not necessarily in phase at any one point of time, while also being co-dependent and interacting over different timescales. A possible source of the mafic material was suggested at a depth of ∼30 km in the mantle (Stern and Benson, 2011), but it is not clear how widespread the source is. Mafic magma influx from this source can, for example, accommodate extension, and thus suppress tectonic (normal) faulting (e.g., Parsons and Thompson, 1991; Buck, 2004; Buck et al., 2006), while high rates of mafic injection into any mush body may overpressure the silicic system, which can in turn cause inflation-induced faulting (Simakin and Ghassemi, 2010; Ruz Ginouves et al., 2021). Here we build upon the concept of three-way silicic-mafic-tectonic interactions presented by Rowland et al. (2010) for the TVZ using more specific observations from Taupō volcano to consider the changing state of volcano-tectonic interactions crossing the interaction zone between Taupō volcano and the rift. In particular, we consider various lines of evidence from the geological, historical and instrumental records to help explain unrest events and eruptive behavior.

From the geological record it is clear that Taupō has a diverse eruptive history that reveals important features of the silicic magmatic system. Rhyolite eruptions through the Holocene vary in volume over several orders of magnitude with no trend in repose time versus eruptive volume (Wilson, 1993), such that the volume of magma erupted does not reflect the amount of time passed since the last eruption, or the repose period until the next eruptive event. This relationship is shared with other global examples of caldera volcanoes (Wilson et al., 2021). These observations imply that the occurrence of small versus large eruptions from Taupō is not simply determined by the steady accumulation of magma into the system through time (Barker et al., 2021), and thus eruptions are likely triggered by processes external to the underlying magma plumbing system. It is likely that magma-tectonic interaction is the critical factor determining the timing and size of many of the eruptive events at Taupō volcano (e.g., Rowland et al., 2010).

This tectonic fingerprint is expressed in several ways. Firstly, individual eruptions have occurred along vent lineaments oriented parallel to the dominant structural fabric of the rift (Wilson, 1993; Houghton et al., 2010). This locus of rift-aligned vents has been repeatedly reactivated throughout the post-Oruanui eruptive history (<25.5 ka), suggesting a supporting role of rift-aligned structures and the WNW-ESE-trending minimum compressive stress (Barker et al., 2021; Figure 2). Secondly, evidence for fault rupture immediately prior to or during eruptions has been described to the north of Lake Taupō (Villamor et al., 2015; Figure 8) and dike-induced rifting has been inferred during the early stages of past silicic eruptions (Allan et al., 2012). Similar patterns and processes have also been documented for Okataina caldera volcano, ∼80 km NNE of Taupo (Nairn et al., 2005; Berryman et al., 2008; Villamor et al., 2011) (Figure 1).

While it is clear that interactions between tectonic and magmatic processes at Taupō volcano can trigger and/or control eruptive episodes, the magmatic system ultimately must be primed for eruption to occur. At Taupō volcano, magma compositions are almost identical between eruptions spaced a few hundreds years apart and of variable volume, but vary over longer time gaps (800–3,000 years) with step-wise shifts in composition due to changing storage conditions (Sutton et al., 2000; Barker et al., 2015). While the processes that form melt dominant magma bodies in any specific silicic system may include multiple mechanisms (e.g., decompression, buoyancy, second boiling, CO₂ flushing; Bachmann and Huber, 2019), the lack of compositional gradients and consistently crystal-poor nature of magmas from Taupō suggests that eruptible rhyolite melt bodies are transient features (Barker et al., 2016), the presence, longevity and extent of which are dependent on currently ill-defined factors.

Despite these uncertainties, the presence or absence of eruptible silicic melt bodies is critical for dictating the response of the system to mafic recharge and/or tectonic processes over various timescales. Indeed, mafic recharge is recognized as a key driver of volcanism (e.g., Leonard et al., 2002; Woods and Cowan, 2009; Rooyakkers et al., 2018), and could also trigger the formation of a melt-rich body (e.g., Barker et al., 2016). Similarly, we suggest that the same logic applies to tectonic processes, which could control the formation or eruption of melt-rich bodies due to static or dynamic stress changes (Davis et al., 2007; Baer et al., 2008; Allan et al., 2012; Biggs et al., 2013; Hamling and Kilgour, 2020; Seropian et al., 2021), or trigger the onset of melt accumulation through, for example, decompression melting.

In addition to triggering eruptions and/or melting episodes, mafic (and occasionally silicic) injections likely trigger many of the seismic crises observed immediately outside of the Taupō caldera (Smith et al., 2007). Intense ground shaking and fault rupture episodes have occurred at the Taupō caldera on 20–40 years timescales in the last 100–200 years, and are often preceded by centimeter to decimeter-scale surface uplifts (Grindley and Hull, 1986). Although there are few data to interrogate the source of these unrest events (e.g., 1922, 1963–1964, and 1983 swarms; Section 5), the nature of the surface deformation cannot be explained by fault-slip alone (Barker et al., 2021), likely suggesting magmatic involvement. For example, during the 2019 unrest, inflation of the magmatic system was the likely cause of seismicity around its margins (Illsley-Kemp et al., 2021), and the magmatic nature of this unrest episode strongly suggests that deep mafic recharge was involved. In all, mafic intrusions at Taupō volcano exert a fundamental control on the dynamics of the magmatic-tectonic system (Barker et al., 2015; Illsley-Kemp et al., 2021).

Drawing upon the described observations from the geological, historical and instrumental records, we consider the sequence of events that could result in unrest and/or eruption at modern Taupō through changes induced by mafic magma intrusion and/or tectonic faulting (Figure 9). Ascent of mafic magmas from depth in the Taupō region is considered to be a relatively common occurrence (perhaps occurring on decadal timescales) (Ellis et al., in review; Barker et al., 2020, 2021; Illsley-Kemp et al., 2021), providing the heat and material required to sustain the large magmatic system (Barker et al., 2021). Tectonic cycles on faults in the region likely occur over longer timescales (10²–10³ years), but may be brought forward (or backwards) by magmatic and/or volcanic processes (Rowland et al., 2010) (see also Villamor et al., 2011, for examples from the Okataina region). In Figure 9, and discussed below, we consider potential event pathways at Taupō volcano with four primary triggers related to short-term changes associated with crustal mafic intrusion and/or tectonic cycling at Taupō:

FIGURE 9. Potential pathways and outcomes of volcano-tectonic interactions at Taupō volcano. Although unrest, and particularly eruptions, can be associated with silicic magma intrusion, these events are ultimately triggered by mafic magma recharge either within the active silicic system (A), outside the area of the active silicic system (B), changes to the tectonic stress state surrounding the volcano (B), or other magmatic processes such as second boiling (A). (A) Within the silicic system, the presence of an eruptible melt body (i.e., >25% melt present; Vigneresse et al., 1996) is key in determining the short-term outcome, but a silicic eruption (Outcome 1) is only possible if the stress state is favorable (e.g., Nostro et al., 1998; Rowland et al., 2010). Eruptible silicic melt body formation, in turn, may be caused by mafic magma recharge, or possibly other internal magmatic processes. Without the presence of an eruptible melt body (i.e., <25% melt present; Vigneresse et al., 1996), the intrusion is likely to stall in the crust due to the presence of the large silicic system and the most likely outcome is unrest (e.g., ground deformation, seismicity: Outcome 2). (B) Outside the active silicic system, mafic intrusions can stall in the crust (Outcome 3), with or without associated surface faulting, or erupt at the surface (Outcome 4). Dynamic and/or static stress changes associated with mafic diking, dike-induced faulting, and tectonic faulting can be imposed on the active silicic system, and can lead to either silicic eruption (Outcome 4) or unrest (Outcome 5) with associated faulting. The short-term outcome is dependent on the presence of an eruptible melt body which, in turn, may be triggered by a change in stress state over longer time intervals. In some cases mafic intrusions away from the silicic system can make their way to the surface to erupt (Outcome 6). Key examples of processes and outcomes are given by the following numbered references: (1) Caricchi et al. (2021); (2) Blake et al. (1992); (3) Allan et al. (2017); (4) Rooyakkers et al. (2018); (5) Barker et al. (2016); (6) Illsley-Kemp et al. (2021); (7) Benson et al. (2021); (8) Rowland et al. (2010); (9) Seropian et al. (2021); (10) Barker et al. (2021); (11) Rowe et al. (2021). See text for further details and discussion.

8.2.1 Mafic Intrusion Within the Modern Silicic Mush System

Mafic recharge may perturb the silicic magmatic system depending on the depth of stalling and the condition of the silicic mush at the time of intrusion. The presence of an eruptible melt body is critical for short-term changes as mafic recharge could push the system to erupt (e.g., Kaharoa, Waimahia; Blake et al., 1992; Leonard et al., 2002). However, we highlight that the system may only erupt if the volcanoes stress field is also favorable for silicic dike propagation to the surface (e.g., Nostro et al., 1998; Rowland et al., 2010; Outcome 1; Figure 9). Over extended timescales mafic recharge may also result in melt body formation or changes in the mush that could eventually contribute to an eruption once triggered. Melt body formation could additionally be formed without the need for mafic intrusion through other internal magmatic processes such as decompression, buoyancy, second boiling, or CO₂ flushing (Caricchi et al., 2021). However, with the frequency of magma-driven unrest at Taupō being 1–2 orders of magnitude more frequent than eruption (eruption frequency of approximately every 500 years during Holocene: Wilson, 1993), it is demonstrably more probable that this pathway leads to unrest (e.g., 2019 seismic swarms, Outcome 2; Figure 9) than any other outcome.

8.2.2 Mafic Intrusion Away From the Silicic Mush System that Stalls in the Crust

Although mafic magmas are inferred to maintain the main silicic reservoir at Taupō, mafic magmas also provide the high heat flow to the other parts of the central TVZ (Bibby et al., 1995), and likely accommodate a significant proportion of crustal extension at depth (i.e., Rowland et al., 2010; Gómez-Vasconcelos et al., 2017). Thus mafic intrusions are likely to be a common occurrence throughout the TVZ, almost all of which stall in the crust with minimal consequence and may escape notice, even in the instrumented age (Outcome 3: Benson et al., 2021; Ellis et al., in review; Figure 9). However, dike-induced stress changes can encourage fault slip (Rubin and Pollard, 1988), and are inferred to have triggered past eruptions from Taupō and other silicic systems globally (e.g., 2005 Dabbahu rifting episode (Ethiopia), 1350 CE Mono-Inyo eruption (United States), 25.5 ka Oruanui eruption (New Zealand); Bursik et al., 2003; Rowland et al., 2007; Allan et al., 2012; Wright et al., 2012). Again, the short-term condition of the silicic system is key when considering the outcome of such a scenario, as a pre-existing melt body could be triggered into eruption by a sudden change in static stress associated with faulting and/or diking (Biggs et al., 2013; Outcome 4; Figure 9). However, the most likely pathway is that the silicic system is not primed and therefore only seismic and geodetic unrest occurs, with the events of 1922 and 1964–1965 being likely examples of this pathway at Taupō volcano (Outcome 5; Figure 9).

8.2.3 Mafic Intrusion Away From the Silicic Mush System that Erupts

The deposits of mafic eruptions are sparse and volumetrically minor (when compared to the rhyolites) in the Taupō area. However, several examples in the geological record exist to the north of the lake (e.g., K-Trig; Figure 2) and are aligned along fault zones in a NNE-SSW orientation (Brown et al., 1994; Leonard et al., 2010). Such eruptions could take place without inducing major faulting and have little effect on the silicic magmatic system (Outcome 6; Figure 9). However, the injection of mafic magma associated with eruptions along these regions could lead to major slip and faulting (e.g., Nobile et al., 2012) that potentially could impact the stability of the silicic magma system to trigger eruption (Outcome 4; Figure 9). Whether mafic intrusion causes large scale fault slip depends on many factors including the location, depth and orientation of the intrusion and the stress state of the surrounding faults (Green et al., 2015).

8.2.4 Tectonic Trigger

In turn, fault slip could simply occur due to normal tectonic loading cycles within the TVZ. Major fault slip within or around the caldera (or even within tens of kilometers) could have significant impacts on the silicic system through changes in dynamic and/or static stresses (Nostro et al., 1998; Manga and Brodsky, 2006; Biggs et al., 2013; Hamling and Kilgour, 2020). Depending on the length extent and geometry of the fault rupture and its proximity to the Taupō magma system, induced static stress changes could create conditions either favorable or unfavorable for magma ascent to the surface (see Nostro et al., 1998), though dynamic stress changes will generally have a destabilizing effect on the magma system (Manga and Brodsky, 2006). The state of the silicic mush system is also a critical factor here, with the presence or absence of an eruptible melt body being a primary control on the outcome. Without a pre-existing eruptible melt body it is most likely that the response of the silicic magma system only leads to unrest (Outcome 5; Figure 9), although changes in stress could assist in the formation of a melt body that may erupt at a later event (Wilson et al., 2021). By contrast, the presence of an eruptible melt body combined with sudden stress changes associated with diking and/or tectonic faulting may produce favorable conditions for generating silicic diking and a possible eruption (Outcome 4; Figure 9). It is also possible that tectonic earthquakes could trigger mafic intrusions which could, in turn, interact with the silicic system (e.g., Chesley et al., 2012).

8.3 Are Volcanic Eruptions or Earthquakes the Primary Hazards at Rift Calderas?

Hosting the most recent supereruption (>450 km³ DRE of erupted material; Miller and Wark, 2008) on Earth (25.5 ka Oruanui eruption) and producing at least 28 further silicic eruptions since, Taupō represents arguably one of the most active silicic volcanoes on Earth. In addition to volcanic hazards, the record of ground shaking and ground deformation in the geological, historical and instrumental records reveals the persistent and ongoing volcano-tectonic hazards posed by this system, both syn- and inter-eruptive. While major syn-eruptive faulting and ground-shaking has been shown to occur, clearly during these time-periods the eruption itself generates the majority of the hazards. Here we consider the hazard posed by volcano-tectonic processes that occur without eruptions.

The 35-year catalogue of seismicity around the volcano indicates an approximately decadal-scale periodicity of unrest at Taupō volcano, with swarms of elevated seismicity occurring every 5–10 years since 1985 (Barker et al., 2021). Some seismic swarms appear to be purely tectonic in origin, but many likely result from the movement of fluids and/or magma (e.g., Illsley-Kemp et al., 2021). Similarly, lake leveling data show that the region outside the lake area is characterized by ongoing vertical surface displacements over 100–200 km² areas, with Fault Belt stations near the northern rim of the lake switching between phases of uplift and subsidence on 10–20 years timescales (Figure 5). Imposed on this near-continuous record of seismicity and surface deformation are periods of intense faulting (i.e., periods of elevated rates of seismicity such as the 1922 and 1983 swarms; Grindley and Hull, 1986; see also Section 5) and associated ground shaking, occurring approximately every 20–40 years (e.g., 1890, 1922, 1964, 1983) and occurring within regions that extend up to ∼15 km from the caldera edge. These events, which can be associated with hundreds of <M3 earthquakes per day (Hull and Grindley, 1984), have produced surface ruptures with discrete centimeter- to meter-scale displacements along multiple fault segments within the Taupō Fault Belt, centimeter-to meter-scale vertical surface movements over broad areas (e.g., 1922 Whakaipo Bay and Kinloch Bay subsidence; see Figure 6 and Section 5.1), ground shaking, and rock sliding.

The overall ∼10 mm/yr extensional strain accumulated in the region also necessitates frequent episodes of seismic energy release. Major faults along the Taupō Fault Belt (e.g., Ngangiho, Kaiapo, Waihī and Whakaipo faults) have time-averaged slip rates on the order of 1 mm/yr and expected recurrence intervals of less than 2000 years (Villamor et al., 2001; Villamor et al., 2015; Stirling et al., 2012), with earthquake magnitudes possibly exceeding M6 for these ∼10 to 20 km-long normal faults (Wells and Coppersmith, 1994; Villamor and Berryman, 2001). An example of this is the >M6 earthquake in 1895, which was potentially sourced from the Kaingaroa Fault 20 km northeast of Taupō volcano (Villamor et al., 2001) (Figure 1).

Silicic eruptions from Taupō have an average recurrence interval of ∼900 years since the 25.5 ka Oruanui eruption, but show no temporal trend with size versus repose period (Wilson, 1993). Seismic crises (e.g., events in 1895, 1922, 1964–1965, 1983, and 2019; see also Section 5) appear to be at least an order of magnitude more frequent. While a silicic volcanic eruption would likely have a greater impact on the central North Island, and possibly the country (Barker et al., 2019), the frequency with which hazardous seismic unrest events are experienced at Taupō highlight earthquake-induced ground shaking as the primary natural hazard to the local populace surrounding the volcano for inter-eruptive periods, i.e., the majority of the time.

Similar to Taupō volcano, earthquakes likely represent the most frequent hazard at other silicic caldera systems situated within regions of crustal extension (or transtension) globally (Gudmundsson, 2020; Acocella, 2021). For example, the Main Ethiopian Rift in East Africa is a region of stretching and thinning of continental lithosphere (∼4 mm/yr extension) and associated magmatism (Casey et al., 2006; Corti, 2009; Rooney et al., 2011; Saria et al., 2014; Birhanu et al., 2016). Silicic centers, situated along-strike every ∼20–50 km, exhibit evidence for Holocene eruptive activity and observed surface deformation suggests the presence of ongoing magma recharge into some of the magmatic systems underlying these volcanoes (Biggs et al., 2011). However, much like Taupō volcano, the frequency with which these centers erupt is at least an order of magnitude less than the frequency of damaging earthquakes. For example, Corbetti volcano exhibits an average eruption recurrence interval of ∼900 years (Martin-Jones et al., 2017), yet the Hawassa rift basin that hosts this volcano has experienced three significant earthquakes (two M >5) since 1960, each associated with rock sliding, and building damage and collapse (Wilks et al., 2017). Similarly, at Long Valley Caldera (United States) persistent earthquake swarms have been recorded within the caldera since the late 1970s (Julian and Sipkin, 1985) and are often attributed to the transport of magma and/or hydrothermal fluids (Julian, 1983; Hill et al., 2003; Shelly et al., 2016). Greater than M5 earthquakes are also documented immediately outside and adjacent to the caldera boundary, with these often exhibiting a foreshock-aftershock sequence (Hill et al., 2003).

Perhaps unsurprisingly the public’s perception of hazards associated with large silicic calderas tends to focus on the catastrophic eruptions of the past. However, our synthesis of volcano-tectonic interactions at Taupō highlights that the primary hazard posed by this volcano, on short-term timescales, is ground-shaking. This trait is likely shared by other silicic calderas in continental rifts worldwide and should be an important consideration for public educators and policy makers in the future. Furthermore, the swarm-like nature of seismicity at Taupō volcano (Grindley and Hull, 1986; Illsley-Kemp et al., 2021; Figure 5) could result in ongoing disruptions to the local populace during unrest events for months to potentially years.

9 Conclusion

We present an overview of data and observations covering a wide range of timescales of deformational processes (1–10⁴ years) at Taupō volcano in order to illustrate how a large silicic caldera volcano interacts with a continental rift. From these observations we show that:

1. Short-term deformation is characterized by decadal-scale uplift and subsidence events with accompanying seismic swarms, ground-shaking and surface ruptures. These events likely reflect recharge into, and shallow crustal intrusion in the vicinity of, the shallow magma mush system (situated at 5–8 km depth).

2. The decadal-scale frequency with which intense seismic events occur highlights the fact that earthquakes, rather than volcanic eruptions, are the primary hazard locally in the Lake Taupō region.

3. Deformation trends on 1–10⁴ yr timescales in the vicinity of Taupō volcano (i.e., within 10 km of the Taupō caldera) are not reflective of the longer-term behavior typical of a continental rift, with crustal magma influx suppressing axial subsidence of the rift basin.

4. Some silicic volcanic eruptions at Taupō volcano are characterized by intense syn-eruptive deformation at least as far as 50 km outside of the caldera structure, including ground shaking, fissuring and fault movements.

5. Eruption and unrest scenarios at Lake Taupō volcano are dependent on the three-way coupling between the mafic-silicic-tectonic system. Unrest events can follow six possible outcomes likely triggered by either mafic injection into or outside the 5–8 km deep magma mush system, or by changes to the tectonic stress state.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://zenodo.org/ doi:10.5281/zenodo.5777477 and fault transect data (10.5281/zenodo.5885833) presented in this study can be downloaded from https://zenodo.org/.

Author Contributions

JM, FI-K, and SB equally contributed to the initial concept and drafting of the paper, based on discussions during eight seminars attended by the coauthors. Figures were constructed by JM, FI-K, SB, EM, RM, MS, PV, PO. PV provided supplemental fault scarp data and analysis of fault extension rates. FI-K, PO, EM compiled and analyzed lake leveling data. JM, CW, FI-K, EM, PV, RM, MS, YS performed supplemental field investigations with logistics support from BS. All authors contributed to manuscript revisions, read, and approved the submitted version.

Funding

This work was conducted and written as part of the ECLIPSE project (contract RTVU1704) funded and supported by the Ministry of Business, Innovation and Employment (New Zealand).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

Melissa Rotella is acknowledged for support organizing and recording these events. Additionally, Jamie Delano, Tim Stahl, and Cecile Massiot attended some events and are thanked for their involvement in these. We are grateful to Gaby Gomez for kindly providing the image for Figure 8C.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/feart.2022.835841/full#supplementary-material

References