Earthquake Analysis Suggests Dyke Intrusion in 2019 Near Tarawera Volcano, New Zealand

Tarawera volcano (New Zealand) is volumetrically dominated by rhyolitic lavas and pyroclastic deposits, but the most recent event in AD 1886 was a basaltic Plinian fissure eruption. In March 2019 a swarm of at least 64 earthquakes occurred to the NE of Tarawera volcano, as recorded by the New Zealand Geohazard Monitoring Network (GeoNet). We use seismological analysis to show that this swarm was most likely caused by a dyke that intruded into the brittle crust between depths of 8–10 km and propagated toward Tarawera volcano for 2 km at a rate of 0.3–0.6 m s−1. We infer that this was a dyke of basaltic composition that was stress-guided toward Tarawera volcano by the topographic load of the volcanic edifice. Dyke intrusions of this nature are most likely a common occurrence but a similar process may have occurred during the 1886 eruption with a dyke sourced from some lateral distance away from the volcano. The 2019 intrusion was not detected by InSAR geodesy and we use synthetic models to show that geodetic monitoring could only detect a ≥6 m wide dyke at these depths. Improvements to geodetic monitoring, combined with detailed seismological analysis, could better detect future magmatic intrusions in the region and serve to help assess ongoing changes in the magmatic system and the associated possibilities of a volcanic event.


INTRODUCTION
Forecasting volcanic eruptions is inherently challenging due to the wide range of unrest signals (e.g., ground deformation, elevated seismicity, gas emissions) that can occur at variable rates and over variable timescales prior to eruption (e.g., Sparks et al., 2012). Many of these unrest signals can also occur without leading to eruption, highlighting the complex nature of the subsurface plumbing systems and the numerous processes that occur beneath dormant volcanoes (Moran et al., 2011). Monitoring subsurface processes is important for developing an understanding of how different volcanoes operate and for assessing any deviation from the normal background state (e.g., Cashman and Sparks, 2013;Acocella, 2014).
Many of the challenges with assessing volcanic unrest are exemplified with caldera volcanoes, which are among the most complex and dangerous types of volcano (Acocella et al., 2015). Caldera volcanoes often cover a wide geographical area (up to several ten of kilometers wide) with large, geometrically complex and heterogeneous magmatic systems that are capable of producing explosive eruptions, sometimes of great size. However, understanding of the current state of these systems is often made complex by the presence of active, shallower hydrothermal systems and fault structures (e.g., Sandri et al., 2017;Mantiloni et al., 2020). Disentangling the signals of magmatic unrest (e.g., variations in the pressure of the system due to magma recharge or ascent) vs. non-magmatic unrest (e.g., changes in the hydrothermal system or tectonic earthquakes) can be difficult, with the potential for many different mechanisms that can produce similar unrest symptoms (Acocella et al., 2015).
Here we document and interpret a series of earthquakes that occurred over a 3 day period near Tarawera volcano (New Zealand) in March 2019. Tarawera is the site of a large basaltic fissure eruption in AD 1886, but is positioned within a larger caldera structure that has a history of producing high-silica (rhyolitic) magmas ( Figure 1). As such, it is important to constrain the origin of these earthquakes and to consider how they relate to the modern magmatic system and ongoing activity in the broader region around the caldera.

GEOLOGICAL BACKGROUND
The Taup o Volcanic Zone (TVZ) in the central North Island of New Zealand (Figure 1) has been volcanically active since approximately 2 Ma, and from 1.85 Ma eruptions have been dominated by high-silica (rhyolitic) magmas (Eastwood et al., 2013;Chambefort et al., 2014). The TVZ can be subdivided into three segments along its length. The southern and northern segments are characterized by andesite volcanism building composite cones, whereas the central TVZ is dominated by voluminous rhyolite volcanism associated with caldera volcanoes and unusually high surface heat-flow (Bibby et al., 1995;Wilson et al., 1995). Over the last 50-60 kyr this rhyolitic volcanism has mainly been focused at Taup o and Okataina volcanoes, at the southern and northern limits of the central TVZ respectively. The TVZ also hosts the Taup o continental rift with present-day rates of extension increasing from ≤5 mm/yr at Ruapehu in the south, to 13-19 mm/yr at the Bay of Plenty coastline and increasing farther offshore (Wallace et al., 2004;Lamarche et al., 2006). The orientation of present-day extension within the Taup o rift varies from rift-orthogonal in the south to oblique in the north (Rowland and Sibson, 2001;Acocella et al., 2003;Townend et al., 2012;Seebeck et al., 2014;Illsley-Kemp et al., 2019). Several studies have suggested that extension within the Taup o rift is partially accommodated by magmatic intrusions (Rowland et al., 2010), however this may not be ubiquitous throughout the rift (Villamor et al., 2011).
The Okataina Volcanic Center (OVC) is built on a system of nested calderas formed by multiple large-volume rhyolite eruptions over at least the last ∼340 kyr (Nairn, 2002;Cole et al., 2010Cole et al., , 2014 (Figure 1). Over the past ∼25 kyr volcanism has been focused at the Haroharo and Tarawera volcanic complexes, in the northern and southern parts of the OVC, respectively (Nairn, 2002). Both complexes are volumetrically dominated by rhyolitic eruptive products, including voluminous lavas forming the respective edifices, erupted from vents aligned along two linear vent zones. A notable feature is that rhyolitic eruptions from the Tarawera complex have often been associated with basaltic intrusion and eruption. The most recent rhyolite eruption (the AD 1314 ± 10 Kaharoa eruption from the Tarawera complex: Hogg et al., 2003) was primed and triggered by the intrusion of primitive basalt (containing primitive olivines: Barker et al., 2020) that was mixed into the rhyolite magma reservoir (Leonard et al., 2002;Nairn et al., 2004). There is also evidence that many of the previous rhyolite eruptions from the Tarawera complex had basaltic input/interaction (Nairn, 1992;Darragh et al., 2006;Shane et al., 2007;Shane et al., 2008). The composition of basaltic magmas erupted from the OVC show little variation between eruptions but they are distinct in their composition compared to basalts erupted outside the caldera to the south (Barker et al., 2020).
The influence of basaltic magmatism at the OVC was most dramatically illustrated in the AD 1886 Tarawera basaltic Plinian eruption, sourced from a ∼17 km long fissure (Nairn, 1979;Walker et al., 1984;Sable et al., 2006) ( Figure 1). However, the 1886 eruptive products show no physical evidence for any interaction with a melt-dominant rhyolite magma reservoir, only shallow xenolithic incorporation (Cole, 1970;Cole et al., 2014;Carey et al., 2007;Carey and Houghton, 2010), although some crustal contamination is indicated from trace element and isotopic data (Gamble et al., 1993;Waight et al., 2017). The inference, therefore, is that the feeder dyke rose from the base of the quartzofeldspathic crust and ascended rapidly while avoiding any interaction with evolved melt-rich magmatic reservoirs. The location of this primitive source, and the pathway of the eruptive dyke are largely unknown. Nairn and Cole (1981) document the 1886 eruptive fissures and find that the surface dykes are arrayed in a left-stepping en-echelon pattern within a larger fissure structure. The dykes have a common orientation of ∼70°, while the eruptive fissure strikes at ∼60°. This led Nairn and Cole (1981) to suggest that orientation of the dyke intrusions reflect the modern day orientation of maximum horizontal compression (SH max ), whereas the eruptive fissure reflects the orientation of an older fault structure. Vents during the eruption began at two locations at the Tarawera summit and spread both NE (for a short distance) and SW (for ∼11 km) along the fissure line during the course of the eruption (Sable et al., 2006).
In addition to the ample evidence that dyke intrusions, specifically of basaltic composition, play an important role in OVC eruptions, they are also thought to play a key role in the accommodation of extension in the Taup o rift (Rowland et al., 2010;Villamor et al., 2011). However, while dyke-accommodated rifting has been directly observed in other continental rifts (Kendall et al., 2005;Ebinger et al., 2013) and realistically has to have occurred during the 1886 eruption, it has not yet been detected in the TVZ in the modern instrumental era. Analysis of GPS data from the OVC area found that strain in the region surrounding the Tarawera fissure should promote the intrusion of dykes (Holden et al., 2015). Other geodetic studies, however, found that the OVC and TVZ are more accurately described as contractional areas at present Dimitrova et al., 2016;Haines and Wallace, 2020).
Knowledge of the modern-day location or state of the magmatic system(s) at OVC is limited. Electrical resistivity imaging shows evidence for a broadly distributed lowresistivity structure across the OVC between 10 and 20 km depth (Heise et al., 2010) and a region of partial melt to the SW of the OVC at depths ≥8 km (Heise et al., 2016). In contrast, seismic anisotropy suggests the presence of a large upper-crustal magma body (Illsley-Kemp et al., 2019). Earthquake activity in the OVC is highly swarm-like, and in the instrumental era (post-1985) has been largely limited to the Haroharo complex and to the southwest of Tarawera beneath Lake Rotomahana, outside the OVC caldera (Hurst et al., 2008;Bannister et al., 2016) (Figure 1). The Lake Rotomahana earthquake activity is thought to be caused by geothermal fluid migration along sub-surface faults between 4 and 7 km depth (Bannister et al., 2016). Eruptive histories and phenocryst mineralogies for the ≤25 ka products of the OVC suggest that Haroharo and Tarawera have different upper crustal magmatic systems, albeit at comparable depths (Cole et al., 2014). Petrological studies based on volatile concentrations in quartzhosted melt inclusions have suggested storage pressures for the erupted rhyolitic magmas in the range of 130-160 MPa (5-7 km) for most examples (Johnson et al., 2011) with greater depths (to ∼12 km) proposed for some (Shane et al., 2007). Also, the presence of cummingtonite as a phenocryst phase in many of the young Okataina rhyolites indicates similar maximum pressures and depths (Nicholls et al., 1992). An improved understanding of the modern system would assist with the monitoring of the volcanoes, particularly as studies by Sherburn and Nairn (2004) and Holden et al. (2017) suggest that it may be possible for significant subsurface magmatic activity to occur in the OVC without detection by the geodetic monitoring currently in place. Building a better understanding of modern magmatic activity is also important for constraining tipping points that could result in magmatic activity/unrest at the OVC cascading toward major unrest and/or eruption (Acocella et al., 2015;Wilson, 2017). In this regard, it is significant that on March 12, 2019, the New Zealand Geohazard Monitoring Network (GeoNet) recorded a cluster of 64 earthquakes to the northeast of Tarawera volcano. In this paper we investigate this seismic cluster and its cause in more detail.

DATA AND METHODS
For this study, we downloaded continuous seismic data from March 12, 2019 to March 15, 2019 (inclusive) from seven shortperiod GeoNet seismometers ( Figure 2). We manually detected earthquakes within the continuous data and picked all visible P Frontiers in Earth Science | www.frontiersin.org January 2021 | Volume 8 | Article 606992 and S phase-arrivals resulting in a total of 348 earthquakes. These earthquakes were initially located with NonLinLoc (Lomax et al., 2000), using a velocity model derived from the nearby Kawerau geothermal field (Clarke et al., 2009) (Supplementary Figure S1). We then used waveform correlation to generate differential pick times with the Obspy package (Beyreuther et al., 2010). We used a 2 s window around each pick, beginning 0.5 s before the pick, allowing the pick to adjust by up to 0.3 s, and each event pair had a maximum hypocentral separation of 8 km. These differential pick times were then used to relocate the entire catalog using the double-difference relocation program GrowClust (Trugman and Shearer, 2017), requiring a minimum correlation of 0.5. This resulted in 94 relocated earthquakes in total ( Figure 2). We calculate local-magnitudes by measuring the peak displacement on a simulated Wood-Anderson seismometer on each seismometer, for each earthquake. We then use the local-magnitude scale developed for New Zealand (Ristau et al., 2016). This results in magnitudes over the range −1.12 to 2.54 M L , and for earthquakes also detected by GeoNet our magnitudes are comparable. We then compute moment tensors using P-wave polarities picked on the wider North Island GeoNet network and the Bayesian moment tensor source inversion software MTfit (Pugh and White, 2018).

RESULTS
The main concentration of earthquakes (131 events total, 71 relocated) occurred in the Puhipuhi embayment to the northeast of Tarawera ( Figure 3) and ranged from −0.47 to 2.54 M L in magnitude. The mean horizontal and vertical errors from the initial location are ±0.8 and ±0.68 km respectively, and the relative horizontal and vertical errors after relocation are ±0.036 and ±0.005 km, respectively. The relocated earthquakes form a clear lineation with a strike of ∼70°. When viewed in cross-section the earthquakes appear to have occurred along a ∼3 km long, vertical plane between 8 and 10.5 km depth ( Figure 3). The earthquakes migrated from NE to SW at an apparent rate of 0.3-0.6 m s −1 , and their range of depths increased at the same time ( Figure 4). These spatial patterns are also produced when using locations generated using GrowClust's internal bootstrapping approach (Efron and Tibshirani, 1991;Efron and Tibshirani, 1994;Trugman and Shearer, 2017), demonstrating that they are not an artifact of data or station distribution (Supplementary Figures S1, S2).
We are able to calculate four reliable moment tensor solutions for earthquakes in this swarm ( Figure 5). These all produce double-couple, normal fault solutions with nodal plane strikes very similar to the strike of the earthquake cluster lineation. The  Figures S4-S7).

Origin of the 2019 Earthquake Swarm
The March 2019 earthquake cluster has several notable characteristics. First, it delineates a vertical structure between ∼8 and 10.5 km depth (Figure 3), which is close to the base of the generally observed seismogenic depth in the TVZ (Bryan et al., 1999). Second, the vertical structure has a strike of ∼70° (  Figure 3), near identical to the observed surficial orientation of dyke intrusions during the 1886 Tarawera eruption (Nairn and Cole, 1981). Third, the earthquakes migrated laterally, from ENE to WSW, at a rate of 0.3-0.6 m s −1 and their depth distribution range expanded during this lateral migration ( Figure 4). There are three possible causes for earthquake swarms in the TVZ; slip on a fault-plane, diffusion of geothermal fluids, or intrusion of a magmatic body. Given that the earthquake cluster occurs below the brittle-ductile transition in the TVZ (Bryan et al., 1999) and occurs along a vertical plane, inconsistent with extensional faulting, we suggest it is unlikely that the cluster was caused by slip along a fault-plane. Lateral propagation of seismicity along vertically inclined faults has been observed at oceanic transform faults and interpreted as aseismic creep along a fault-plane (Roland and McGuire, 2009). While these observations are similar to those presented here, the tectonic setting in the TVZ is quite different and we feel it is therefore unlikely that our observations are explained by strike-slip creep. Smith et al. (2007) proposed that subsidence at Taup o caldera (New Zealand) in 1983 was caused by the dewatering of a magma body at depth, and at Yellowstone caldera (USA) this process has been shown to cause propagating seismicity on a vertical plane, similar to our observations (Waite and Smith, 2002). However, the diffusion velocity of hydrothermal fluids resulting from dewatering is  Figure 3. This shows that the earthquakes increase in their depth extent through time. Bottom: Earthquake distance along section A (Figure 3) vs. time. This shows a clear propagation of the seismicity WSW along-strike at a rate of 0.3-0.6 m s −1 .
Frontiers in Earth Science | www.frontiersin.org January 2021 | Volume 8 | Article 606992 5 estimated, and observed, to be several orders of magnitude slower than the intrusion rate we observe, and thus our data is not satisfactorily explained by the diffusion of geothermal fluids from a dewatering magma body. In addition, upwelling geothermal fluids have been observed to cause seismicity in nearby Rotomahana (Bannister et al. 2016; Figure 1). However, the seismicity at Rotomahana shows no lateral migration, occurs between 4 and 6.5 km depth, and occurs in repetitive "bursts" over at least the past 25 years (Bannister et al., 2016). These characteristics suggest that the March 2019 earthquake cluster is unlikely to be caused by a similar upwelling of geothermal fluids.
We therefore suggest that the 2019 earthquake swarm was most likely caused by a dyke, intruding to depths as shallow as ∼8-10.5 km, that migrated toward Tarawera volcano in a WSW direction. The migrating earthquakes are inferred to have been caused by country rock opening at the tip of the propagating dyke (Belachew et al., 2011;Sigmundsson et al., 2015), and this is evidenced by the normal fault moment tensors ( Figure 5). These moment tensors may also be explained by a component of non-double-couple deformation (Supplementary Figures S4-S7), which would also support a dyke model. We can thus use the migration rate of earthquakes as a proxy for the intrusion rate of the dyke (i.e. 0.3-0.6 m s −1 ). Our calculated intrusion rate of 0.3-0.6 m s −1 is very similar to those recorded for basaltic dyke intrusions in the Afar rift (Belachew et al., 2011;Barnie et al., 2015), Red Sea (Eyles et al., 2018), El Hierro (Martí et al., 2013), and Iceland (Einarsson and Brandsdóttir, 1980). However, it is significantly faster than intrusion rates calculated from a rhyolite dyke intrusion in the Main Ethiopian Rift (3 × 10 −3 m s −1 ; Temtime et al., 2020). We therefore propose that the dyke intrusion was most likely basaltic in composition.

Volume Estimate
Dyke intrusions represent movement of magma that can often be detected by volcanic monitoring through measuring surface deformation. However, the ability to detect intrusions depends on many factors including the volume and location of the intrusion. Based on the energy balance caused by an intrusion and the seismic moment release, Bonaccorso et al. (2017) derived an empirical relationship relating seismic moment to dyke thickness. Using this relationship, and the total seismic moment release of 1.77×10 13 N m, the thickness of the 2019 dyke near Tarawera would be ∼0.15 cm. Alternatively, global observations of dyke intrusions suggest length to thickness ratios of 1,500 (Gudmundsson, 1987(Gudmundsson, , 2011 giving a thickness of ∼2 m. Therefore, based on the dyke's length and associated seismic moment release, we could expect an intrusion with a volume of 1.1-15 × 10 6 m 3 . Ground deformation as a result of discrete dyking events has been well documented in a number of volcanic and rift environments (Hamling et al., 2009;Sigmundsson et al., 2015;Hamling et al., 2019). To investigate whether the inferred dyke was detected geodetically, we use ascending and descending synthetic aperture radar (SAR) data acquired by the Sentinel-1 satellite to generate two co-intrusion interferograms ( Figure 6). Neither of the interferograms show an obvious signal consistent with an intrusion, and there is no detectable signal in the GNSS network ( Figure 2). However, given the depth of the dyke intrusion, it is possible that ground deformation was too small to be detected.
To investigate how much volume could be hidden before deformation would be detected by the InSAR data, we generated a suite of synthetic intrusions. We fixed the geometry of the dyke based on the earthquake observations and examined the effect of the dyke thickness on the detectable signal. For each thickness value, we generated a synthetic interferogram by adding the modeled deformation to the ascending and descending interferograms respectively. We then evaluated the signal to noise ratio using the deformation from the modeled intrusion as the reference signal. For both the ascending and descending datasets, the signal to noise ratio only becomes positive once the thickness of the dyke exceeds ∼5-6 m ( Figure 6). This suggests that for a similar sized intrusion, at these depths, we would be unlikely to detect it with current geodetic monitoring unless its volume exceeds ∼37-45 × 10 6 m 3 .
Our estimated intrusion volume of 1.1-15 × 10 6 m 3 , based on dyke length and seismic moment release, should be considered as a minimum estimate, as the maximum depth of the FIGURE 5 | Moment tensors for four earthquakes associated with the dyke intrusion, earthquake ID numbers are below each. Locations are indicated by white stars in Figure 3. Events 2019p191131, 2019p191139, and 2019p191176 occur within the deep dyke-associated cluster, while event 2019p191204 occurs in the shallow cluster at a depth of 1.3 km. Moment tensors are calculated using MTfit (Pugh and White, 2018) and show the full range of possible nodal planes, with the highest probability solution shown in green. Further details of the moment tensor solutions is shown in Supplementary Figures S4-S7. Circles and triangles represent positive and negative polarity picks respectively.
Frontiers in Earth Science | www.frontiersin.org January 2021 | Volume 8 | Article 606992 6 earthquake sequence may not provide a lower limit on the intrusion. Evidence from elsewhere in the TVZ (e.g., Bryan et al., 1999), suggests the depths of the earthquake sequence may only reflect the maximum depth to which brittle failure occurred. It is therefore likely that the source reservoir for the 2019 dyke intrusion lies at depths ≥10 km in the vicinity of the Puhipuhi embayment. Further geophysical imaging in this region may be able to determine the location and depth of this magma reservoir.

Structural and Magmatic Implications
Our results lend contemporary evidence to the suggestion of Villamor et al. (2011) that extension within the OVC is accommodated by dyke intrusions. We also observe shallow (1-4 km depth) earthquakes directly above the dyke during intrusion (Figure 3), this is most likely caused by dykeinduced mechanical faulting accommodating extension in the brittle shallow crust (e.g., Belachew et al., 2011). Dyke intrusions can also yield good indications of the orientation of the crustal stress field (Wadge et al., 2016). The orientation (∼70°) of this dyke intrusion suggests that this is also the orientation of the maximum horizontal compression (SH max ). This orientation differs slightly to observations from the rest of the TVZ which suggest an SH max orientation of 40°-60° (Townend et al., 2012;Illsley-Kemp et al., 2019). It has been previously suggested that the direction of extension in the Taup o rift becomes more oblique in the Okataina region (Seebeck et al., 2014), possibly influenced by the presence of a large magma reservoir (Ellis et al., 2014). Our results support these previous observations and may suggest a higher degree of rift-obliquity within the Okataina region than that proposed by Seebeck et al. (2014). However, multiple studies have also shown that dyke orientations and trajectories can be influenced by local changes to crustal stresses. Topographic loads (e.g., volcanic edifices) have been shown to modify the local stress field and cause dykes to propagate toward them (Maccaferri et al., 2011;Rivalta et al., 2015), whereas topographic depressions (e.g., calderas) promote circumferential dyke intrusions (Gaete et al., 2019). The 2019 dyke intrusion propagated toward the Tarawera volcanic edifice, perpendicular to the major OVC caldera boundary, but sub-parallel to the Puhipuhi embayment margin (Nairn, 2002, Figure 3). For comparison, the 3.4 ka Rotokawau basalt in the northwestern region of the OVC (Figure 1) erupted along a fissure that is perpendicular to both the OVC caldera and Haroharo volcanic complex, and oblique to the rift orientation (Nairn, 2002). Shortly following the Kaharoa eruption there were multiple hydrothermal eruptions SW of Tarawera. Nairn et al. (2005) propose that these were caused by the intrusion of a SW-NE oriented basaltic dyke. Similarly, the 1886 eruptive fissure extended along-axis to the SW beyond the OVC FIGURE 6 | Line of site displacement for ascending (top) and descending (bottom) SAR data acquired by Sentinel-1. The data (left) shows no observable ground deformation from the dyke intrusion. We then use synthetic dykes, centered on the 2019 dyke location but with varying thickness, to show that only a dyke intrusion with thickness ≥6 m would be detectable.
Frontiers in Earth Science | www.frontiersin.org January 2021 | Volume 8 | Article 606992 caldera margin and into Lake Rotomahana, showing no evidence that its extent or trajectory were influenced by the OVC caldera boundary (Nairn, 1979;Nairn, 2002). This suggests that within the OVC, the load of the volcanic edifices (rather than the collapse caldera margin) are exerting a first order control on crustal stress and dyke trajectories, promoting dyke propagation toward the base of these volcanic edifices (Maccaferri et al., 2011;Rivalta et al., 2015). These are important considerations for volcanic hazard assessment and for interpreting future signals of unrest in the OVC. Given that dyke intrusions are shown to propagate toward the volcanic edifices from outside the OVC caldera, we speculate that the basalts which have been shown to influence Tarawera eruptions are not necessarily sourced from directly beneath the volcano. A lateral migration of basaltic magma toward Tarawera may also help explain the lack of evidence for interaction with the silicic magma reservoir in the 1886 eruption (Carey et al., 2007;Carey and Houghton, 2010). In such a case, the 1886 dyke may have completely bypassed the silicic system that is presumably still present at depth (but not yet detected by geophysics) beneath Tarawera from the Kaharoa eruption in AD 1314. Alternatively, the 1886 dyke may have intersected the silicic system beneath Tarawera, but the silicic system had since cooled beyond the solidus and is now moribund. Further geophysical assessment of the magmatic plumbing system beneath the OVC is required to address these unknowns.
It is likely that the 2019 dyke intrusion may have had similarities to the early dyke intrusions that fed the 1886 Tarawera eruption. But what were the controlling factors that caused the 2019 dyke to arrest, whereas the 1886 dyke intrusion led to a deadly and violent eruption? There are numerous potential causes for the arrest of the 2019 dyke, including magma solidification within the dyke (Fialko and Rubin, 1999), pressure decrease at the source magma reservoir (Rivalta, 2010), or stress and/or structural barriers. It is also important to consider that although the earthquake activity ceased, the 2019 dyke may have continued to propagate aseismically due to the crust becoming more ductile in the region closer to the caldera boundary. If a dyke intrusion were to reach the Tarawera Volcanic complex, we suggest that the load of the Tarawera edifice would cause it to shallow (Rivalta et al., 2015) and it could then take advantage of pre-existing crustal weaknesses to propagate to the surface and erupt. This may have been the case in 1886, and if a dyke intrusion were to propagate toward, and reach, Tarawera volcano in the modern day then a similar basaltic eruption may occur. However, it's important to note that geological evidence suggests that a high flux of primitive melt is required to sustain the rhyolite volcanism observed in the TVZ (Barker et al., 2020). This means that dyke intrusions are likely to be common and thus the vast majority do not lead to an eruption. Developing a better understanding of the tipping points that cause dyke intrusions to erupt on rare occasions is an important research goal.

Implications for Monitoring and Future Activity
Our analysis of a series of earthquakes that occurred in 2019 suggests a recent dyke intrusion that propagated toward Tarawera volcano, the site of New Zealand's largest and most destructive historical eruption (Walker et al., 1984). The OVC is monitored by GeoNet, primarily through seismometers and GNSS sensors though they are fairly sparse in the Puhipuhi embayment ( Figure 2). In an active magmatic system, dyke intrusion can occur as a common process without necessarily resulting in an eruption (Acocella, 2014). However, detecting these events is critical for evaluating volcanic hazards and for monitoring ongoing surface changes. If a future dyke intrusion was identified early in its propagation, and the dyke dimensions were constrained by geodetic techniques, then the expected total mechanical energy release could be calculated using derived relationships (Bonaccorso et al., 2017;Bonaccorso and Giampiccolo, 2020). This, alongside detailed seismic analysis of the dyke-induced earthquakes, would allow for the forecast of the time and location of dyke arrest, and the associated likelihood of significant unrest and/or eruption (e.g., Aspinall et al., 2006;Constantinescu et al., 2016). We have shown that the requisite seismological analysis is possible with the current GeoNet seismic monitoring system albeit with detailed analysis long after the earthquake swarm. This process could be significantly improved using real-time detection and location of earthquakes (e.g., Chamberlain et al., 2020). However the geodetic monitoring network, particularly GNSS station density, would need to be improved in order to detect and constrain a similar future dyke intrusion. We have also shown that dyke intrusions sourced from outside the OVC caldera can propagate toward Tarawera, and we speculate that similar "external" dyke intrusions may have influenced past eruptive activity. Therefore, monitoring of Tarawera volcano should not focus solely on the region directly beneath the volcanic edifice but consider the broader surrounding area.

CONCLUSION
On the 12th March 2019, a swarm of 131 earthquakes occurred to the NE of Tarawera volcano, within the OVC. We use highprecision earthquake locations to demonstrate that the swarm was most likely caused by the WSW intrusion of a ∼2 km long dyke, with seismicity occurring between 8 and 11 km depth. By tracking the migration of seismicity we estimate that the dyke propagated at a rate of 0.3-0.6 m s −1 , a very similar propagation speed to that observed for basaltic dyke intrusions in global extensional tectonic settings. The intrusion rate, coupled with the depth of the dyke intrusion, indicate that this was most likely a basaltic dyke intrusion. Based on the dyke length and seismic moment release we estimate a total (minimum) intrusion volume of 1.1-15×10 6 m 3 , and use a suite of synthetic models to show that the associated ground deformation could not have been detected by co-intrusion interferograms. The 2019 dyke intrusion propagated toward the Tarawera volcanic complex from outside the OVC, suggesting that this topographic load has a first-order control over the local crustal stress. Improvements to geodetic monitoring in the OVC could help to better identify future dyke intrusions and allow for the early assessment of the volcanic hazard they may or may not pose.
Frontiers in Earth Science | www.frontiersin.org January 2021 | Volume 8 | Article 606992 8 DATA AVAILABILITY STATEMENT All seismic data used in this study are freely available from GeoNet (https://www.geonet.org.nz/). The earthquake catalogue is provided in QuakeML format in the Supplementary Material and as an online dataset: 10.5281/zenodo.4035171.

AUTHOR CONTRIBUTIONS
TB performed seismological analysis. FI-K supervised TB and led writing of the article. HE produced Figure 1 and helped write the article. IH performed geodetic analysis and helped write the article. MS supervised TB and helped write the article. CW helped write the article. EM supervised TB and helped write the article. SB helped write the article.

FUNDING
TB was funded by a Victoria University of Wellington Summer Scholarship and the ECLIPSE program, which is funded by the New Zealand Ministry of Business, Innovation and Employment. FI-K, EM are fully funded, and SB, CW, IH, and MS are part supported by the ECLIPSE program. HE is part supported by the ECLIPSE program and a Victoria University of Wellington doctoral scholarship.

ACKNOWLEDGMENTS
This project analyzed and plotted the seismic data using EQcorrscan (Chamberlain et al., 2018), Obspy (Beyreuther et al., 2010), GMT (Wessel et al., 2019), and Pyrocko (Heimann et al., 2017). We thank Nico Fournier for constructive feedback on the manuscript. We also thank the constructive reviews of Yohei Yukutake and Francesca Di Luccio.