https://orcid.org/0000-0001-5753-6568
SUMMARY
Magnetotelluric data from Mount Tongariro have been analysed using an unstructured tetrahedral finite-element inversion code that incorporates topography, which was not included in previous analysis of these data. Incorporating topography adds information, which stabilizes the resistivity inversion modelling, and for the first time allows details of the shallow hydrothermal system and its relationship with the underlying magmatic system to be resolved. Specifically, an electrically conductive zone between 4 and 12.5 km depth marks the underlying magmatic system, which is shown to directly connect via conductive pathways to the area where the most recent phreatic eruptions at Tongariro occurred in 2012. The resultant phreatic eruptions in 2012 August and November showed no new magmatic component to the eruption deposits. Nevertheless, by combining the magnetotelluric resistivity image with relocated seismicity, we can see that seismicity (a proxy for magma ascent) migrated from the top of the magmatic system into the hydrothermal system in the months preceding these eruptions. Magmatic interaction with the extant hydrothermal system likely caused the over pressurization for the phreatic eruption. This work highlights the utility of combining geophysical methods and petrological data to constrain phreatic eruption processes.
INTRODUCTION
The Tongariro Volcanic Centre in the central North Island of New Zealand is a large andesite–dacite stratovolcano complex located at the southern termination of the Taupo Volcanic Zone (TVZ). The TVZ is a rifted arc formed by the subduction of the Pacific plate along the east coast of the North Island (Fig. 1) and is characterized by a very high geothermal heat flux (Bibby et al. 1995) and high rates of magma generation; ∼10 000 km3 of predominately rhyolitic magma having been erupted from the TVZ over the last ∼2 Ma (Wilson et al. 1995).
Inset map shows the tectonic setting of New Zealand. Red lines mark the position of the Australia-Pacific plate boundary and Hikurangi subduction margin (red triangles). The orange box in the insert shows the location of the shaded relief map of Mount Tongariro that forms the rest of the figure. Yellow triangles show the locations of the three most recently active eruption vents (Ngāuruhoe, Red Crater and Te Maari). Lavas < 25 ka old are indicated in pink (after Pure et al. 2020). MT measurement sites are shown by red and green dots; and the green dots show the location of sites for which the model fit is shown in Fig. 6. Coordinates of the southeastern corner are given in New Zealand transverse Mercator projection.
At the southern end of the TVZ, where rifting and volcanism cease, seismic tomography indicates that the brittle crust thickens from 25 to 40 km (Reyners et al. 2006) and there is an absence of mantle return flow. In contrast, return flow has been interpreted to be the underlying cause of the central TVZ's high heat-flux and the high-rate of magmatism (Reyners et al. 2006, Eberhard Phillips et al. 2020).
The southern termination of the TVZ is marked by three active andesitic volcanoes: Mt Ruapehu, Mt Ngāuruhoe and Mt Tongariro, collectively referred to as the Tongariro Volcanic Centre. The oldest lavas from Mt Tongariro are dated to be ∼350 ka old (Pure et al. 2020) while dates from Ruapehu are younger (∼200 ka; Gamble et al. 2003).
Mounts Tongariro and Ngāuruhoe, the focus of this study, form a single composite stratocone volcano of dacitic–andesitic composition (e.g. Leonard et al. 2021). At least 15 eruption vents have been recognized as being active on Mt Tongariro in the last 250 000 yr, all of which lie within an extensional graben (Villamor et al. 2017; Milicich et al. 2020) aligned SW-NE from Tama lakes to Te Maari (Fig. 1).
An episode of large Plinian to sub-Plinian eruptions from multiple vents occurred during a period of about 200 yr between 11.2 and 11.4 ka. These eruptions, known as the Mangamate episode (Nairn et al. 1998), were composed of mainly hydrous hornblende–dacite tephras (Auer et al. 2016; Heinrich et al. 2020, 2022). Although the mafic melt source for these eruptions is inferred to reside in the upper mantle (at >20 km depth), a more silicic (i.e. intermediate composition) mush zone between 5 and 12 km depth is inferred to be the source of the Mangamate eruptive products (Auer et al. 2016; Heinrich et al. 2022).
After this unusual eruptive episode, the eruptions at Tongariro returned to their pre-Mangamate style characterized by numerous small (<0.1 km3), short-lived (≤1 ka) eruptions composed of two-pyroxene plagioclase andesites thought to originate from depths between 4 and 21 km (Hobden et al. 1999; Moebis et al. 2011; Arpa et al. 2017). Historic Tongariro Volcanic Centre eruptions (i.e. within the last 150 yr) have occurred from Red Crater, Mt Ngāuruhoe and Te Maari (Fig. 1) and most recently from Mt Ngāuruhoe in 1977 and from the Upper Te Maari Craters on the northern side of Mt Tongariro in 2012 August and November (Hurst et al. 2014; Scott & Potter 2014).
The 2012 eruptions are interpreted to have been phreatic and triggered by a dyke intrusion leading to over-pressurizing of the shallow hydrothermal system as the ash samples do not appear to contain juvenile magma (Jolly et al. 2014; Pardo et al. 2014; Miller et al. 2018). The source depth of these small phreatic eruptions is inferred to be shallow, that is, ∼500 m depth (1100 m asl) consistent with the depth of the hydrothermal system (Hamling et al. 2016) although triggered by deeper magmatic activity.
Since an interconnected network of conductive fluid or melt lowers the conductivity of normal crustal rocks by orders of magnitude, magnetotelluric (MT) measurements provide a useful method for imaging zones of magma storage (partial melt) and/or associated magmatic fluids. At shallow levels, zones of clay-rich hydrothermal alteration, geothermal brines and acidic fluid are also conductive and are imaged well by MT. Throughout the central TVZ, extensive MT surveys have shown conductive bodies in the crust that are interpreted as deeper intermediate mush zones or shallow bodies of crystallizing magma and high-temperature (saline) fluid (Heise et al. 2007, 2010; Bertrand et al. 2012, 2013, 2015, 2016, 2022).
At the southern termination of the TVZ, an early MT study of Mt Ruapehu (Ingham et al. 2009) identified a dyke-like feature, that provides a pathway for gas and fluid beneath the crater. However, the connection to a deeper magmatic source could not be determined by these data, due to their limited spatial coverage. At Mt Tongariro, an extensive 3-D MT survey indicated crustal magma at ∼4–12 km depth (Hill et al. 2015), with measurements repeated at 14 locations following the 2012 Te Maari eruptions (Hill et al. 2020).
A limitation of the 3-D inversion modelling technique used in previous studies, is that topography could not be included. In areas of pronounced topography, such as Mounts Ruapehu and Tongariro, short-period MT data will be strongly affected by topographic effects (Jiracek 1990), compromising resolution of near surface structure, and potentially introducing artefacts in the inversion models. Uncertainty in the depth of deeper resistivity structure introduced by the assumption that the earth is flat also limits effective comparison with the distribution of seismicity. By adequately including topography in the modelling procedure, additional prior information is added which helps constrain and stabilize the inversion.
Here we use a 3-D finite-element inversion code FEMTIC (Usui 2015; Usui et al. 2017) that includes topography to re-analyse MT data acquired in 2009 and 2010 at Mount Tongariro (Hill et al. 2015). This inversion code represents the topography and subsurface using an unstructured tetrahedral mesh. The model mesh used here was designed to optimize the resolution of both the shallow hydrothermal system, and the depth range suggested for magma storage. Outside the area of data coverage, the size of the mesh elements is coarsened to reduce the computational burden. In this paper, we combine our resistivity inversion model of Tongariro MT data using FEMTIC with accurately re-located seismicity to image magma ascent prior to the most recent phreatic eruption at Te Maari. We also explore different approaches to galvanic distortion analysis introduced into MT data by near-surface resistivity heterogeneities.
MT DATA
The MT method measures the natural variations of the earth magnetic field and the induced electric field over a broad period range. If the MT measurements are distributed over a wide area, they allow the subsurface electrical conductivity (or resistivity) structure to be inferred down to depths comparable to about half the aperture of data coverage, provided that the period range is also sufficient. The depth to which the electromagnetic signals penetrate depends on the bulk electrical resistivity (ρ) of the region and the period (T) of the MT signal. For a uniform half-space, the skin depth (δ) is given by δ ≈ 0.5 |$\surd $|ρT in km, where T is in seconds and ρ is in Ω·m. This quantity can be used (with caution) as a general indicator of the detection depth in more complicated situations.
A dense array of 128 broad-band MT measurements with a site spacing of ∼1 km were measured on Mt Tongariro (locations shown in Fig. 1) between 2009 and 2010; before the 2012 eruptions at Te Maari. These data were all recorded using Phoenix Geophysics MTU-5 and 5A data loggers and induction coils. Electric field measurements were made using 70-m grounded dipoles and Pb-PbCl2 electrodes. At most locations reliable data were recorded in the period range between about 3 ms and ∼1000 s. The MT survey covered an area of 600 km2, that is, ∼30 km wide in the east–west direction and ∼20 km in the north–south direction, allowing for good resolution to depths of ∼15 km.
By normalizing the measured electric field with the magnetic field in the frequency domain, MT data can be converted into a set of period-dependent transfer functions that express the magnitude and phase relationships between the electromagnetic field components. These transfer functions are normally expressed by a 2-D impedance tensor (units Ω) and an induction vector that describes the relationship between the horizontal and vertical magnetic field components.
The phase relationships inherent in the impedance tensor data are usually re-expressed as the MT phase tensor (Φ), which can be represented graphically as an ellipse (Caldwell et al. 2004). Importantly, maps showing phase tensor ellipses and induction dvectors allow the main structural features present in MT data to be discerned before inversion and provides a convenient method of data visualization that is largely free of the distorting effects (galvanic-distortion) cause by near-surface conductivity heterogeneities. For example, the long axis of a phase tensor ellipse at a measurement site located on the conductive side of a 2-D resistivity structure is oriented parallel to the strike of the structure. In contrast, the ellipse is oriented perpendicular to the structure if the measurement site is located on the resistive side. In a 3-D setting, the ellipses tend to take on parallel or perpendicular orientation to the strongest resistivity gradient within a skin depth of the measurement period.
In Fig. 2, we show phase tensor ellipse maps for periods of 0.05, 0.3, 2.7, 21 and 170 s, with the ellipse colour indicating the value of the geometric mean (Φ2) of the principal axes (Φmin and Φmax). Values of Φ2 (< 45 °) shown by cold colours indicate increasing resistivity with depth, whereas high values of Φ2 (> 45°), indicated by warm colours show increasing conductivity. In these maps, the ellipse sizes have been normalized by Φmax. At short periods (Figs 2a and b), high Φ2 values mark the line of eruptive vents and show the conductive alteration associated with the hydrothermal system. At longer periods (Figs 2c and d) ellipses located east of the eruptive vents are strongly polarized and align parallel to the vent axis with high Φ2 values, indicating a very strong resistivity contrast with respect to the western side, which must be more resistive. At the longest period (Fig. 2e), the phase tensor data indicate increasing resistivity structure at depth, but still show a distinct difference east and west of the vent axis, consistent with the eastern side being more conductive.
(a)–(e) MT data shown as phase tensor ellipses and induction vectors at periods of 0.05, 0.3, 2.7, 21 and 170 s. The colour fill shows Φ2, with values > 45° indicating increasing conductivity with depth with values < 45° indicating decreasing conductivity. Arrows show induction vectors using the sign convention that vectors point towards regions of greater conductance. Black lines on the eastern side of the data coverage show high voltage AC transmission lines. Shaded relief map as in Fig. 1.
At periods between 0.05 and 0.3 s, the induction vectors on both the eastern and western sides point at the line of vents indicating the presence of a narrow conductive zone. At increasing periods (2.7 and 21 s), a large contrast between the eastern and western half of the survey area is observed, with the induction vectors on the western side pointing south-west towards the relatively conductive sediments in the Whanganui basin (e.g. Naish et al. 1998). At longer periods (>170 s, Fig. 2e), the induction vectors in the east have small magnitudes indicating the presence of a deeper conductor. Note that the short-period induction vector data at the easternmost sites are not shown in Figs 2(a)–(c) as these data were strongly affected by high-voltage power transmission lines (e.g. Bertrand et al. 2012) that run just to the east of these measurement locations. Induction vector data at periods < 20 s for these locations were also not included in the subsequent inversion modelling.
3-D RESISTIVITY INVERSION MODELLING
3-D tetrahedral mesh
Creating a finite-element mesh that is fine enough to both accurately represent topography and resolve features at the depths of interest, without being too computationally large is a crucial first step for inversion success. To find an optimum mesh, we used a topographic grid with 200 m sampling and a commercial meshing package called Cubit (Coreform Cubit 2021). Beyond a 24 × 34 km rectangle centred at the median of the MT measurement locations, the topographic grid was smoothed and tapered to a constant height at a distance of 100 km from the centre of the grid.
For our 30 × 20 km MT data aperture with site spacing (∼ 1 km) and period range 0.0125–341 s, the range where the data sensitivity is good, is between approximately 0.25 and 20 km depth. The skin depth of the data at long periods will be greater than 20 km but without lateral control from a wider area, the structure at depth is less well constrained.
For the inversion modelling, we used a starting model with 100 Ω·m resistivity and air resistivity of 109 Ω·m. At the longest period included (340 s), for 100 Ω·m the skin depth is 92 km less than the distance to the coast of the North Island. For this reason, the coastline and oceans surrounding the North Island were not included in the inversion modelling, simplifying the mesh design.
The total size of the mesh volume was 800 km × 800 km laterally, and 400 km thick, with the upper 100 km representing air above the topography. Our mesh requires a large variation in cell size, which we accomplished by making extensive use of Cubit's Exodus II-based field function. At the ground surface in the vicinity of the observation locations, we specify a background resolution of 1 km. We then define a sphere of 600 m radius surrounding each site, with a specified size of 60 m at the centre of the sphere and a size of 150 m at the outer surface of the sphere. We also specify half ellipsoids above and below the ground surface, with target cell sizes of 1 km within each ellipsoid. Outside of the ellipsoids, the target cell sizes increase rapidly to a value of 80 km near the edges of the mesh. Details of the mesh are shown in Fig. 3. The depth radius of the subsurface ellipsoid was determined by trial and error to give an element edge length of about 1 km at about 12 km depth, resulting in a mesh with ∼700 000 elements (Fig. 3).
(a) Triangular surface facets of the tetrahedral mesh elements forming topographic surface. The dashed line shows the location of the NW-SE cross-section through the finite-element mesh shown in (b). Surface facets for the entire mesh are shown in (c).
Mesh generation was an iterative process, where we first started with a mesh having relatively uniform element sizes. We then created sizing information for this mesh and performed the meshing process again to yield a mesh with a better element size distribution. We finally performed the same process a third time using the second iteration of the mesh, yielding the mesh used for this study. Each iteration of mesh generation took approximately 1 hr on a standard Linux workstation.
Inversions using a finer (3000 000 elements) and coarser (450 000 elements) mesh were also tested. Although these meshes resulted in similar resistivity models, the coarse mesh did not provide adequate spatial resolution in the depth range (between 5 and 10 km) where magma storage was expected to occur. Compute times using the fine mesh made regularization parameter testing and trials of different approaches to galvanic distortion (discussed later) burdensome, with insignificant gain in final resolution. The results shown in this paper were computed using four 512 GB memory nodes, 4 tasks per node and 16 Xeon Gold 6338 cores @2.00 GHz per task on GNS Science's hpc cluster. Each iteration took approximately 15 min to complete.
Input data
Input data for the inversions comprise the components of the impedance tensors and induction vectors along with their corresponding errors at 16 periods between 0.0125 and 341 s. Data quality is generally good (e.g. Fig. 2). Isolated data outliers were also removed before inversion. In an area northeast of Mount Ngāuruhoe, phase tensor maximums are large compared with the minimum phase (see Fig. 2d) In this region, the phase tensor is badly conditioned and sensitive to small magnetic-sensor alignment errors and data noise, which can result in negative phase tensor determinants. Data points with negative determinants were also omitted. At periods less than 20 s, noise generated by high-voltage AC transmission lines can strongly affect nearby vertical magnetic field measurements (e.g. Bertrand et al. 2012). For this reason, only the induction vector data with periods >20 s have been included from the 14 measurement sites located within ∼5 km of the transmission lines on the eastern side of the survey area (Fig. 2).
Error floors
The error floor of all impedance tensor elements was set at 3 per cent of the geometric mean (|$\sqrt {|Zxy} Zyx|$|) of the off-diagonal impedance-tensor components. We also tested an impedance error floor of 5 per cent for the inversion which resulted in a normalized rms value < 1, suggesting that the inversion was under constrained. The error floor used for induction vector data was set at 0.015 in both cases.
Model regularization
FEMTIC uses a single regularization parameter (α2) to control the trade-off between model roughness and the data fit measured by the normalized root mean squared (rms) misfit. Larger values of α2 result in smoother models but also larger rms misfit values. To determine the inversion model that best balances the competing demands of obtaining a smooth model with an acceptable data fit, inversions with α2 values of 2, 4, 8, 16, 32, 64 and 128 were calculated. A plot of the resulting L-shaped trade-off curve is shown in Fig. 4. These models converged (i.e. the rms misfit ceased to change significantly) after 21, 11, 10, 9, 7, 6 and 6 iterations, respectively. In the discussion of different approaches to the galvanic distortion that follows, the inversion models were all computed using an α2 value of 8.
L-shaped trade-off curve for different regularization parameter values (α2 = 2, 4, 8, 16, 32, 64 and 128). The black dot marks the position of our preferred model (α2 = 8, overall normalized rms = 1.22) on the L-curve.
Galvanic distortion
Galvanic distortion of the measured MT data is caused by near-surface inhomogeneities that are too small to be resolved by the measurement spacing. Bibby et al. (2005) describe a method to calculate the form of the distortion, represented by a real valued 2 × 2 matrix (a tensor), based on a prior phase tensor analysis of the impedance data. Where this analysis shows that a section of the impedance data at short periods has a 1-D behaviour, the distortion matrix can be calculated to within a single multiplicative scale factor known as the static shift. In this circumstance, the distortion can be stripped out of the measured impedance and the (non-unique) scale factor chosen to be the magnitude of the measured (perhaps shifted) 1-D impedance at short periods.
Inverting impedance data where distortion has been removed adds the prior information that the resistivity structure of near-surface elements around the stripped sites is 1-D. However, if the static distortion is caused by the topography rather than a localized heterogeneity near the measurement point, this will be inconsistent with the topographic response at short periods. The inversion will then produce near-surface structure around the measurement point to compensate. Other approaches to the distortion problem are possible, most obviously by including the elements of the distortion matrix at each measurement site as unknowns to be recovered by inversion (e.g. Avdeeva et al. 2015); an option also implemented in FEMTIC.
In Fig. 5, we show phase tensor pseudo-sections for the uniform resistivity (100 Ω·m) initial model used for the inversions. In the area of highest topography, Φ2 values (Fig. 5a) deviate significantly from the uniform resistivity half-space value of 45°. In this area, the topographic effect is significant, and the phase perturbation persists to periods of about 10 s. However, at most sites the difference between phase tensor principal values (Fig. 5b) is small (< 2°) at periods greater than about 0.1 s. This means that at most sites the phase response produced by the topography is 1-D to a good approximation and the effect of a uniform resistivity topography at long periods can be represented as a galvanic distortion and static shift.
Phase tensor pseudo-section for the uniform 100 Ω·m resistivity topographic model for the locations shown by green dots in Fig. 1. Colour filling of the tensor ellipses and contours in (a) show Φ2. The difference between the phase tensor principal axes (Φmax − Φmin) are shown in (b).
The finite-element mesh of the topography was constructed using a 100 m × 100 m digital terrain map. In areas of very steep terrain, this mesh will only partially capture the topographic effect at measurement sites located adjacent to abrupt changes in height (e.g. sites located near the edge of a crater wall). At these sites the residual galvanic distortion caused by the topography will be stripped out of the data if the method of Bibby et al. (2005) is applied.
Three inversion models using different approaches to distortion analysis were tested. First, we inverted the stripped impedances by removing galvanic distortion using the method described in Bibby et al. (2005), which was possible at 93 of the 128 sites available. This model is shown in Figs 6(a) and (d); with the inversion algorithm converging after 10 iterations with a normalized rms misfit of 1.23. Second, we included the elements of the distortion matrix as additional unknowns for the inversion to calculate (Figs 6b and e). Inverting for the distortion matrix increases the number of unknowns and converged more slowly reaching a final rms misfit of 1.25 after 14 iterations. Finally, we inverted the original (unstripped) impedances. This inversion model converged after 11 iterations with an overall rms misfit of 1.3 (Figs 6c and f).
Surface element resistivities for impedance and induction vector data inversions using different approaches to galvanic distortion analysis. Contours of the measurement site rms misfit are shown in white. (a) Inversion using stripped impedance tensor data (see the text). (b) Inversion including distortion tensor elements at each site as additional unknowns to be determined. (c) Inversion using the original (unstripped) impedance data. The magenta line shows the location for the cross-sections. (d) Cross-section with TNG076 modelled with the (e) stripped data, calculated data and (f) original (unstripped data). (g) Site TNG076, lines show the stripped data. Map limits in (a)–(c) are the same as in Fig. 1.
At deeper levels, the differences between the inversion models are indiscernible but, as is shown in Figs 6(d)–(f), the near-surface resistivity structure is noticeably different near some of the measurement sites. Comparing data misfit from each measurement site (contours in Fig. 6) shows that the largest differences occurred at only a few points that have large static shifts and lower data quality (Fig. 6g).
The overall normalized rms misfit is greatest for the inversion model that used the measured data and lowest for the stripped data model, although the difference is small. This small difference suggests that fine meshing adjacent to measurement points provides the additional degrees of freedom needed to adequately model the distortion without prior analysis. That said, while the difference is minor, the model that inverted the stripped data has the smoothest near surface resistivity structure and the smallest rms misfit, so is our preferred model.
To illustrate the quality of the data fit in our preferred inversion model, Fig. 7 shows the observed and modelled phases tensors for 13 sites located along an NW-SE profile (location shown in Fig. 1). Greatest misfits (contours in Fig. 7d) occur in the area NE of Mount Ngāuruhoe where the phase tensor ellipses for periods ∼20 s are elongated, the phase-tensor skew angles (not shown) are large (>30°), and the data are most sensitive to noise and magnetic sensor misalignment errors. Models with less smoothing (α2 = 4 and 2) reduce the misfit in this area, but these models are significantly rougher. However, it is clear by comparison of the measured and calculated phase tensors in Fig. 7 that there is a very good overall fit to the measured data along this profile, through the centre of the MT array.
(a) Phase tensor ellipses at all periods used in the inversion at the locations shown by green dots in Fig. 1, ellipse colour showing Φ2. (b) Calculated phase tensors from the preferred model. Contours show the rms 0.5*√(Δxx2+ Δxy2+ Δyx2+ Δyy2) of the tensor misfit (Δ) defined as Δ = I—Φcal *Φ−1obs, where I is the identity matrix and Φobs and Φcal are the observed and calculated phase tensors.
Seismicity
Seismic events for this study were identified using the deep-learning-based EQTransformer algorithm, which detects and picks P- and S-phase arrivals utilizing a pre-trained model (Mousavi et al. 2020). We applied EQTransformer to 12 yr (2010–2021) of continuous streams of three-component data recorded by all permanent GeoNet seismometers (Hurst et al. 2014) in the vicinity of Mt Tongariro and Ruapehu. The P and S arrivals identified using EQTransformer were then used to derive initial hypocentre locations using NonLinLoc (Lomax et al. 2000), which locates individual events using a probabilistic, nonlinear global search. Final event hypocentres for 1435 events were then derived using the double-difference tomoDD algorithm (Zhang & Thurber 2003), for all events with more than five phase observations, using event-pair phase-time differential times together with the absolute phase arrival times to refine relative event locations. Event-pair differential catalogue phase times were calculated from the absolute phase data for all pairs of events separated by less than 10 km. We used a fixed 3-D velocity model interpolated from the regional model of Eberhart-Phillips et al. (2022) for traveltime calculations.
For events shallower than 10 km, the number of locatable events derived using EQTransformer was similar to the number of events in the GeoNet catalogue; 779 events were derived using EQtransformer in the 30-km radius around Tongariro in a 12-month period in 2022–2023, compared with 739 in the GeoNet catalogue for the same time period. No magnitude threshold was applied—events were not located if there were less than 6 (P or S) phase observations, which effectively filters out smaller events that are not well observed.
3-D INVERSION RESULTS—RESISTIVITY STRUCTURE
Fig. 8 shows a map-view slice through our preferred resistivity model at 4 km depth with epicentres of seismicity between 2.5 and 5.5 km depth superimposed. A narrow dyke-like conductive zone (resistivity < 30 Ω·m) is present between Te Maari and Red Crater, aligned with the axis of the eruption vents. Seismicity in this depth range is concentrated at the NE margin of the conductive zone near Te Maari and in a small zone about 2 km to the SW of Te Maari within the conductor.
Map view of a depth slice at 4 km of the 3-D resistivity model. Triangles show the locations of the three most recently active eruption vents. White dots show earthquakes between 2.5 and 5.5 km depth. The locations of the NW-SE and the SW-NE profiles in Fig. 10 are shown by white dashed lines.
Figs 9(a) and (b) show orthogonal NW-SE and NE-SW cross-sections (locations shown in Fig. 3a) through the resistivity model, crossing at Te Maari. Seismicity located within 1.5 km of each cross-section has been superimposed. In Fig. 9(a), a wine-bottle-shaped conductive zone, tilted to the NW, extends upwards from about 15 km depth reaching the surface at Te Maari. Seismicity is concentrated in the neck of this conductive zone at about 4 km depth.
(a) NW-SE cross-section through the Te Maari vent. Earthquakes, within 1.5 km of the cross-section, are shown by white and black dots; black indicating the earthquakes that occurred in the 6-month period before the 2012 August eruption. (b) SW-NE cross-section along the line of vents. Vertical white dashed lines show where the two cross-sections intersect.
The high conductivity near the surface between Red Crater and Te Maari corresponds to areas of clay-rich hydrothermal alteration and acidic fluids. The hydrothermal system extends beneath Mount Ngāuruhoe although the surficial part of the cone appears to be unaltered, consistent with the magnetic anomaly observed there (Miller & Williams-Jones 2016).
As can be seen in Fig. 9(b), and more clearly in Fig. 10, at about 1.5 km beneath the surface, resistivities increase except in the Te Maari area, where the conductive zone extends downward to join the deeper conductive body at about 4.5 km depth. Below Te Maari, seismicity is concentrated in this narrow conductive zone linking the deep conductor to the surface. In particular, seismicity in the 6-month period before the 2012 eruption (shown by black dots in Figs 9a and b is concentrated here and is noticeably shallower than prior seismicity shown by white dots in Figs 9a and b).
3-D perspective views of the most conductive regions of the resistivity model. (a)–(d) Show faceted iso-surfaces enclosing volumes of progressively increasing resistivity. Black dots show the MT measurement site locations at the surface.
DISCUSSION
We interpret the earthquakes as being located above a large body of partial melt, shown by the deeper conductor between 4.5 to 14.5 km below sea level. This conductor is approximately 3 km × 3 km wide and dips to the northeast (Figs 9 and 10). The most conductive part of this body (< 5 Ω·m, Fig. 10) occurs below the Te Maari area, which we interpret as the volume where the melt fraction is greatest.
To estimate the melt fraction, the chemical composition of the melt and the absolute resistivity of the magma body must be known (Gaillard 2004). However, MT best resolves resistivity gradients rather than the resistivity. This factor along with the conductive hydrothermal system at the near surface, decreases our ability to determine the bulk resistivity at depth.
In Fig. 11(f), we show the forward responses for the MT yx-polarization (comprised of measured variations in north–south magnetic fields and east–west electric fields) at three sites for different resistivity models (A–E) of the melt body. These sites are located at the margins of the melt body at depth, where the differences in the response between models are greatest (Figs 11a–e). In contrast, the responses for the xy-polarization (not shown) are indistinguishable.
(a)–(e) 8 km depth slices for the resistivity models in Table 1. (a) Resistivities of all the mesh elements within the 10 Ω·m iso-surface in Fig. 10 set to 1 Ω·m. (b) Element resistivities within the 5 Ω·m iso-surface set to 1 Ω·m. (c) Corresponding depth slice through the preferred inversion model. (d) and (e) Element resistivities within the 10 and 20 Ω·m iso-surfaces set to 10 and 20 Ω·m, respectively. (f) Corresponding yx apparent resistivity and impedance-phase curves for these models (different colours) at the locations shown by the white dots in (a)–(e).
In model A (Fig. 11a), all resistivities within the 10 Ω·m iso-surface in Fig. 10 were changed to 1 Ω·m. This change increases the bulk conductance of the region within the iso-surface compared with the inversion model (Figs 8, 9, 10 and 11c) where the resistivity lies in the range from 3 to 10 Ω·m. This model has the largest misfit and does not fit the data at long periods. In model B (Fig. 10b), all resistivities within the 5 Ω·m iso-surface in Fig. 10 were changed to 1 Ω·m. The bulk conductance of this model is still greater than our preferred inversion model (Fig. 11c), but has a similar misfit and fits the data within the error bars (Fig. 11f). The conductance was then decreased with respect to the inversion model. In model D (Fig. 10d), all resistivities within the 10 Ω·m iso-surface (Fig. 10) were changed to 10 Ω·m and in model E (Fig. 10e) all resistivities within the 20 Ω·m iso-surface (Fig. 10) were changed to 20 Ω·m. Model D almost fits the data within the error bars whereas model E, which has a larger rms, does not.
These results are summarized in Table 1 where the models have been ordered by decreasing bulk conductance from top to bottom. Our preferred inversion model lies in the middle of the conductance range tested and has the smallest rms.
Rms misfit values for test models. Models B and D also fit the data within the error bounds.
Rms misfit values for test models. Models B and D also fit the data within the error bounds.
CALCULATION OF MELT PERCENTAGE
The melt conductivity (i.e. the conductivity of the liquid fraction) is strongly influenced by the content of Na2O and H2O (Gaillard 2004). At Tongariro, the deep source (primitive) magmas are H2O poor, ∼ 1–2 wt. per cent (Lormand et al. 2020), but as they cool and crystallize, the H2O content of the melt will increase (e.g. Blake 1984), decreasing the resistivity. At Tongariro, various authors (Hobden et al. 2001; Auer et al. 2016; Arpa et al. 2017; Pure et al. 2020) have measured H2O and Na2O content of the melt inclusions between 1–7 and 2.3–3.7 per cent, respectively. Calculation using the melt percentage calculator SIGMELTS (Pommier & LeTrong 2011) with these chemistries produces melt resistivities between 0.3 and 10 Ω·m.
To estimate the melt fraction, we used Archie's law (Archie 1942) ρ = a ρf ϕ−m where ρ is the bulk resistivity, ϕ melt fraction and ρf melt resistivity and the empirical parameters suggested by ten Grotenhuis et al. (2005) for low melt fractions (a = 1.47 and m = 1.30). Assuming the lower bound on the melt resistivity (0.3 Ω·m) based on the chemistry and the bulk resistivity (10 Ω·m) suggested by our hypothesis testing, gives a melt fraction of 9 per cent. More resistive melts necessarily result in higher melt fractions (e.g. ϕ =50 per cent for ρf = 3 Ω·m and ρ= 10 Ω·m). Using a different inversion algorithm, Hill et al. (2015) suggest a melt fraction of 18–45 per cent for a magma body of (bulk) resistivity 0.2–0.7 Ω·m. The difference between Hill et al. (2015) and our estimates of the melt fractions highlights the difficulty of melt fraction calculation based on MT inversion results.
The volume of the conductor in Hill et al. (2015) is significantly smaller than the ∼50 km3 conductor in our model (Figs 8 and 9). To recover approximately the same total conductance within a smaller volume, Hill et al. (2015) require much lower resistivities. So, in this (limited) sense, the inversion results are roughly consistent. Also, Hill et al. (2015) used a larger error floor (20 per cent for the off-diagonal elements and 10 per cent for the diagonal elements) and a different regularization in the horizontal and vertical directions. For our results, a single error floor was used (3 per cent of √(|Zxy) Zyx|) for all of the impedance tensor components; a significantly stronger constraint on the data than used by Hill et al. (2015). A further difference between the inversions is that we included the induction vector data, which helped better constrain the deeper structure. However, we believe that the most important factor in terms of the difference between the models is the introduction of topography (prior information is being added) with the least important factor being the difference in approach taken to distortion removal.
Auer et al. (2016) argue that the small magma batches (e.g. Hobden et al. 2002; Lormand et al. 2020) that characterize Tongariro's eruptions before and after the dacitic Mangamate eruption are accompanied by an overall decrease in the water content of the melt during storage and ascent. If so, the most conductive part of the MT model (Fig. 10a) may represent a hydrous intermediate composition magma body that now resides in the mid-crust. We speculate that the presence of a more silicic body, with relatively low density, deflects ascending (possibly more mafic) magma batches to the periphery of the mid crustal conductor (i.e. Maccaferri et al. 2014) as illustrated in Fig. 9.
A striking feature of the mid crustal conductor is that it appears to dip to the northeast away from the vent area to link with a deeper conductor at about 22 km depth. A conductor at this depth is plausible and could be interpreted as mantle melts ponding near the base of the crust. Based on the composition of eruptive products from recent eruptions at Red Crater, Gruender et al. (2023) present evidence of pre-eruptive recharge of a primitive magma from a deep-mantle source, which is consistent with this interpretation. However, this conductor is not well resolved due to the limited aperture of the MT survey (20–30 km).
The shallowing earthquakes before the initial 2012 eruption suggest that magma, gases, and fluids ascended into, and over-pressurized the hydrothermal system (e.g. Hurst et al. 2014; Jolly et al. 2014; Hill et al. 2020). The over-pressured hydrothermal system was subsequently un-roofed after a portion of the overburden failed, generating a polymict debris avalanche (Proctor et al. 2014).
CONCLUSIONS
A 3-D finite-element inversion model that includes topography has allowed the hydrothermal system and its connection to the deeper magmatic system at Mount Tongariro to be better resolved than previous studies. A key factor in the increased resolution was a careful mesh design that included fine structure at the depths of interest, but coarsening outwards to also allowed fast computation time, which facilitated efficient hypothesis testing and trying different approaches to distortion analysis. These analyses confirm that if the mesh adjacent to the measurement sites is fine enough, the inversion will model the distortion without prior analysis. However, by inverting the stripped data we have implicitly added prior information to the inversion, which produces a smoother surface structure and a marginally smaller rms value.
The 3-D inversion model shows a remarkable correlation with the seismicity and a compelling image of the magmatic system broadly consistent with petrological and geological evidence. In particular, the seismicity prior to the 2012 eruptions appears to show the ascent of magma from a deeper chamber into the hydrothermal system, resulting in the 2012 phreatic eruptions at Te Maari. Our results also support the conclusion of Heinrich et al. (2022) that Mt Tongariro is capable of large Plinian style eruptions and that a large silicic magma chamber is present.
ACKNOWLEDGEMENTS
This work was funded by the Geothermal: The Next Generation research programme funded by the New Zealand Ministry of Business, Innovation and Employment (MBIE contract C05X1904). This work was also supported by MBIE through the Hazards and Risk Management and Understanding Zealandia/Te Riu-a-Māui Programmes (Strategic Science Investment Fund, contract C05X1702).
DATA AVAILABILITY
Data sets related to this paper can be requested at https://doi.org/10.21420/0ky8-mx13, hosted at the GNS Science Dataset Catalogue.
