Research ArticleGEOPHYSICS

Off-axis magmatism along a subaerial back-arc rift: Observations from the Taupo Volcanic Zone, New Zealand

See allHide authors and affiliations

Science Advances  03 Jun 2016:
Vol. 2, no. 6, e1600288
DOI: 10.1126/sciadv.1600288


Continental rifting and seafloor spreading play a fundamental role in the generation of new crust. However, the distribution of magma and its relationship with tectonics and volcanism remain poorly understood, particularly in back-arc settings. We show evidence for a large, long-lived, off-axis magmatic intrusion located on the margin of the Taupo Volcanic Zone, New Zealand. Geodetic data acquired since the 1950s show evidence for uplift outside of the region of active extension, consistent with the inflation of a magmatic body at a depth of ~9.5 km. Satellite radar interferometry and Global Positioning System data suggest that there was an increase in the inflation rate from 2003 to 2011, which correlates with intense earthquake activity in the region. Our results suggest that the continued growth of a large magmatic body may represent the birth of a new magma chamber on the margins of a back-arc rift system.

  • InSAR
  • New Zealand
  • Volcanology
  • Stress transfer
  • Taupo Volcanic Zone
  • Off-axis magmatism


Our understanding of rifting episodes and emplacement of magma beneath oceanic spreading centers is constrained by a few, albeit well-documented, events (1, 2). Observations of subaerial rifting episodes in Afar (3, 4) and Iceland (3, 5, 6) have begun to shed light on the complex interactions between volcanism, tectonics, and magmatism in regions of incipient seafloor spreading. However, in arc and back-arc settings, there are very few candidate rift zones that produce significant volumes of magma and are observable above sea level. The subaerial Taupo Volcanic Zone (TVZ), New Zealand, is the most productive area of silicic volcanism on Earth (7) and offers a unique opportunity to investigate the role of magma in the development of complex back-arc rift.

The TVZ is located in the central North Island of New Zealand and formed as the result of continued subduction of the Pacific plate beneath the Australian plate at rates of 38 to 49 mm/year (Fig. 1) (8, 9). Extension across the TVZ over the last 1.6 million years (My) is expressed at the surface by extensive northeast-southwest–trending fault belts along a narrow ~30-km-wide zone (10) with active volcanism during the last 300 thousand years (11). Previous studies have split the TVZ into three distinct segments (southern, central, and northern TVZ) on the basis of variations in volcanism along the rift (11, 12). In the central TVZ, volcanism over the last ~2 My has led to the eruption of more than 10,000 km3 of predominately rhyolitic magma and resulted in the formation of eight calderas (11). Conductive bodies at depths of 6 to 10 km are thought to be zones of interconnected melt (13), consistent with seismic velocity models, which suggest a heavily intruded shallow crust beneath the central TVZ at depth (14). Large-scale subsidence along the length of the central TVZ, observed by interferometric synthetic aperture radar (InSAR) and Global Positioning System (GPS) data, has been interpreted predominantly as the cooling and contraction of magma at a depth of ~6 km (15). There has been no volcanic activity in the central and onshore northern TVZ since the 1886 Tarawera eruption when a ~17-km-long basaltic fissure eruption occurred along the southern margin of the Okataina Caldera (Fig. 1) (16). At the North of the Okataina Caldera (Fig. 1), rifting is considered to be amagmatic with extension being accommodated by normal faulting (12). Across the northern TVZ, swarms of earthquakes have been detected since the development of the seismic network in the late 1970s. However, the nature of the swarms, whether tectonic or magmatic in origin, has remained largely unresolved (17, 18).

Fig. 1 Location maps showing the North Island of New Zealand and the study area at the northern end of the Taupo Volcanic Zone.

(Left) Color-shaded relief map of the North Island of New Zealand indicating the study area (black box) and boundaries of the currently active TVZ [<300 thousand years ago (ka); solid blue line] and old TVZ [>300 ka; dashed blue line). (Right) Color-shaded relief of the study area highlighting the location of the modern TVZ (gray transparency), Okataina and Rotorua Caldera boundaries (blue dashed lines), and the relocated 2005–2009 Matata earthquake swarm, which are color-coded by year. The dashed box indicates the region shown in Fig. 4. The moment tensor (MT) solution for the largest event is shown by the beach ball. The inset in the top left shows the cumulative number of earthquakes during the sequence.


Using synthetic aperture radar (SAR) images acquired by the European Space Agency (ESA) Envisat satellite between 2004 and 2011, we generate an InSAR time series covering the northern TVZ (Figs. 1 and 2, fig. S1, and Methods). Line of sight (LOS) increases of up to 12 mm/year are observed at the south of Rotorua within the central TVZ (Fig. 2 and fig. S2), consistent with previous observations of subsidence explained by the cooling magma at depth (15). We also identify LOS increases of ~15 mm/year centered over the Kawerau geothermal field within the Rangitaiki Plains (Fig. 2 and fig. S2). To the west of the currently active rift, adjacent to the Bay of Plenty coast, we observe up to 10 mm/year of LOS decrease over a ~300-km2 area (Fig. 2). Two continuous GPS (CGPS) sites (Fig. 2 and fig. S3), operating since around 2007 (operated as campaign sites in 2005 and 2007), and one campaign site with observations since 2007 indicate uplift and up to 10 mm/year of southward motion, consistent with the InSAR observations. InSAR and GPS time series show that the deformation has been ongoing since at least 2004, but in late 2010, after the final InSAR acquisition, both of the GPS sites show a decrease in the uplift rate (Fig. 2 and fig. S3).

The magnitude and spatial pattern of deformation along the Bay of Plenty coast are not consistent with a tectonic origin, and given the proximity to the active TVZ, we assume that the observed uplift is a result of the accumulation of magma at depth. Despite the apparent volcanic origin of the deformation, the broad pattern of uplift cannot be explained by the intrusion of a dyke. We therefore model the inflation as a point source (19) embedded in an elastic half-space to represent the accumulation of magma at depth. In addition to the Envisat track and GPS velocities, we use the data from the Advanced Land Observing Satellite (ALOS) track 324 presented by Hamling et al. (15) (Fig. 3 and fig. S4). We solve for the position, depth, and volume change of the source by performing a hybrid Monte-Carlo downhill simplex inversion, which minimizes the square misfit between the observed and calculated displacements at each observation point (20). To assess the distribution of model parameters, we perform a bootstrap resampling procedure to generate 1000 new data sets and reinvert each new data set (fig. S4). The best-fitting model indicates a source located just offshore at a depth of ~12.9 ± 2.1 km inflating by ~0.01 ± 0.006 km3/year. Although the point source approximation provides a good fit to the data, we also model the deformation, assuming a more realistic sheet-like source (sill) (21), which would require a much lower overpressure to inflate (22). For example, inflating a preexisting magma chamber with a radius of 2 km by 0.01 km3 would require a pressure increase of ~12 MPa. For an equivalent sill (10 × 10 km with 10 cm of opening), the required pressure would be ~0.4 MPa. We assume that the sill is horizontal and solve for the body’s position, depth, strike, length, width, and amount of inflation. As with the Mogi source, the distribution of model parameters is determined by the inversion of 1000 bootstrapped data sets (fig. S5). The best-fitting model (Fig. 3 and fig. S5) suggests a 19.1 ± 3.0 × 24.2 ± 3.9–km source, inflating by ~20 ± 13 mm/year at a depth of 9.5 ± 2.1 km, producing an annual volume increase of ~0.009 km3, comparable to the Mogi source. Trade-offs between different model parameters and large uncertainties in the position and size of the body are a function of the data distribution (fig. S5). With most of the observations limited to the Bay of Plenty coast, the offshore position and the extent of the magma body are difficult to assess.

Fig. 2 Surface velocities derived from InSAR, leveling, and GPS data over the Bay of Plenty coast.

(A and B) InSAR time series for scatterers (red circles) colocated with the two CGPS sites (RGOP and RGMT) along with the GPS velocities converted into the satellites LOS. (C) Mean LOS velocities derived using 32 Envisat SAR images as described in the main text. The black lines indicate mapped active faults. The black dashed circle shows subsidence associated with the Kawerau geothermal power plant. The black, white, yellow, and red arrows show the velocities for the two continuous (RGOP and RGMT) and campaign (EAY8) GPS sites (circles and lines show the 95 and 68% confidence levels for horizontal and vertical velocities, respectively) for the two periods described in the main text. The colored triangles show the GPS displacements (15) at campaign and continuous sites converted into the satellite’s LOS. The blue arrows show the velocities derived from leveling data collected in 1950 and 1978. The image shown in the black box shows the errors associated with the LOS displacements shown in the main figure.

Fig. 3 Observed, modeled, and residual velocities based on the best-fitting tensile dislocation model.

The black and red arrows show the observed horizontal and vertical velocities, respectively. The white and yellow arrows show the modeled horizontal and vertical velocities, respectively. The black and white arrows in the residual panels on the right hand side of the figure show the horizontal and vertical residuals, respectively. The black dashed line marks the edge of our best-fitting sill source, and the black dot shows the location of the Mogi source. The red dashed line in the first panel shows the boundary of the present-day TVZ.


The location of the observed uplift lies outside of the present-day boundary of the active TVZ (Fig. 1) (11). Although the region falls within the limits of the older TVZ, evidence for an off-axis source is surprising, considering that the last known volcanism in the area occurred some 425 ka (23). The northern region of the TVZ is thought to be largely amagmatic with the accommodation of long-term extension through normal faulting (12). Numerical simulations of dyke intrusions show that rift topography is an important factor in the transfer of magma from depth (24). For nonvolcanic and wide rifts, models suggest the formation of stacked sills in the lower crust and, in some cases, off-axis volcanism if intrusions continue to propagate to the surface. Other models suggest that off-axis volcanism develops due to the interaction of boundary faults and magma at depth (25). Small displacements along the bounding normal faults cause dilatational stresses to build at the base of the fault, promoting the intrusion of magma. Although there has been no recent volcanism in the area, the location and depth of our best-fit source (Fig. 4 and fig. S6), located on the footwall side of the inferred bounding faults, fit with these models for the formation of an off-axis magma body.

Fig. 4 Stress change and b value analysis of the Matata earthquake swarm.

(A) Coulomb failure stress (CFS) change calculated for each of the earthquake epicenters of the 2005–2009 Matata earthquake swarm based on our best-fit source model (gray rectangle) and rake of −100°. (B) Estimate of the b value for the Matata earthquake swarm as described in the main text based on a node radius of 2 km. (C) Graph showing the percentage of events within the swarm with positive stress changes against rake direction for our best-fit Mogi and sill sources. The histogram shows number of events and their rake (45° bins) for events with MT solutions.

Although the injection of magma into the crust without subsequent eruption is not uncommon along volcanic arcs (26, 27), the development, and eventual eruption, from large sill complexes imaged beneath supervolcanoes is thought to be a result of the incremental accumulation of magma over long periods of time (28). It is possible that the inflation event observed here is not an isolated incident. Geological estimates of uplift based on the elevation of stranded beach ridges suggest rates of ~3 mm/year over the last ~1700 years (29). GPS data in the region, most notably at EAY8, continue to show uplift of 4 to 5 mm/year since the end of 2010 (Fig. 2), with lower rates of uplift observed at RGMT and RGOP (fig. S3). Furthermore, leveling data collected between 1950 and 1978 along the Bay of Plenty coast (30) (Fig. 2) show uplift of ~100 mm (~4 mm/year) across the same region. Assuming the same sill source as described earlier, this would require an average annual inflation rate of ~0.0025 km3/year since the 1950s, 35% of the inferred rate from 2004 to 2011.

Repeated swarms of earthquakes have been detected in the Bay of Plenty region since 1977, including a well-recorded swarm between 2005 and 2009. Increased seismicity associated with the intrusion of magma into the crust has been observed at a number of volcanic centers (31). To investigate whether the timing of the 2005–2009 Matata earthquake swarm and the observed deformation are related, we use our best-fit source model to compute the Coulomb failure stress for each of the earthquakes within the relocated catalog. For each event, we generate a small rectangular fault, whose center corresponds to the epicenter of the event, and fix the strike to the trend of the swarm (38°). Using the formulations of Okada (32), we estimate and reproject the stress at the epicentral location onto the plane for a range of rake directions (Fig. 4 and fig. S6). MTs for some of the larger events in the swarm [weight-average molecular weight (Mw) > 3.4] are predominantly normal with a small component of strike-slip (Fig. 4 and fig. S7) (17). We find that ~99% of the events are associated with increases in the Coulomb stress for rakes of between −70° (pure normal) and −140° (oblique normal), matching the available MT solutions (Fig. 4 and fig. S7).

The b value, which is a measure of the relative number of small to large earthquakes that occur in a given area over a given time, has been shown to vary in both space and time on the basis of the local tectonic conditions (33). High b values (>1) are often associated with volcanic regions due to local heterogeneities, the presence of fluids, and high temperatures (34). To test whether there is any spatial variation in the earthquake behavior, we compute the b value associated with the 2005–2009 Matata earthquake swarm (Methods). We find that the largest b values (>1.5) are located close to the coastline, in the vicinity of our inferred source, and decrease offshore (Figs. 4 and fig. S8). Although we estimate increased stresses for nearly all of the earthquakes in the swarm, the larger b values close to the inferred source suggest that the events early in the sequence were triggered as a result of the intruding magma. These events may then trigger subsequent seismicity further offshore where we see a reduction in the b value.

With only subaerial observations, it is impossible to assess the extent of the deformation offshore. However, two small islands located ~12-km offshore show uplift of ~4 mm/year over the observation period (fig. S2). Marine gravity data (35) indicate a negative anomaly coincident with the location of our best-fitting source, and bright reflectors observed along a seismic line acquired ~25-km offshore have been interpreted as magmatic sills at a depth of ~10 km (36), suggesting that the source may extend offshore. If the source does extend offshore, then our volume estimates provide a lower bound for the amount of magma being accumulated. Between 2004 and 2011, we estimate that ~0.06 km3 of new material was intruded along the Bay of Plenty coast. If the inflation rate from 1950 to 1978 is representative of the time in between, then an additional ~0.2-km3 magma may have been emplaced since 1950. We suggest that long-term uplift along the Bay of Plenty coast is the result of the continuous intrusion of magma at depth and that the occurrence of frequent earthquake swarms in the region is in response to periodic increases in melt supply, like that observed between 2005 and 2011. The city of Tauranga is located only ~50 km to the west of the activity, emphasizing the need to better understand rifting processes in the region for improved volcanic and seismic hazard assessment. Continued monitoring will allow us to identify future pulses of inflation and potentially enable us to predict earthquake swarms in the region. Although the ultimate fate of the magma remains unclear, its presence may represent the birth of a new magma chamber on the margins of arguably the world’s most active region of silicic volcanism, which has witnessed 25 caldera-forming eruptions over the last 1.6 My (7).


To analyze the SAR data, we followed the method of Hooper et al. (37) using StaMPS. Thirty-two Envisat SAR data (fig. S1) were focused using the Jet Propulsion Laboratory/CalTech ROI PAC software (38), and each interferogram was generated using the Doris software (39). Topographic corrections were made using a 3–arc sec (90 min) digital elevation model generated by the NASA Shuttle Radar Topography Mission (40). To reduce the number of noisy points selected, points with SDs of more than 1 radian were removed. Residual orbital errors were estimated and removed by fitting a quadratic plane from each interferogram in the time series. To account for variations in the stratified water vapor content, we estimated and removed the best-fitting stratified delay, assuming a linear function with height. Using the corrected displacement time series, we solved for the best-fitting displacement rate x such thatEmbedded Imagewhere the design matrix A contains the time interval of each epoch relative to the master epoch (6 December 2005), b is a matrix containing the displacements at each scatterer for each epoch, and Σ is the variance-covariance matrix. In total, we identified ~53,000 individual scatterers across the region (Figs. 2 and fig. S2). The GPS data were analyzed using the GAMIT/GLOBK software version 10.6 using over 100 CGPS stations in the analysis to evaluate site positions in the ITRF08 fixed Australian reference frame. In the analysis, we solved for station coordinates, satellite orbit, and Earth rotation parameters, estimating atmospheric zenith delays every 2 hours and using three atmospheric gradients per day. We used the IGS08 azimuth and an elevation-dependent absolute phase center model with an elevation cutoff angle of 10° for the ground-based antennas and applied the FES2004 ocean loading model. Before becoming continuous sites in 2007, RGOP and RGMT operated as campaign sites with measurements in 2005 and 2007, which were used in estimating the velocity for the same observation period as the InSAR data. RGOP was also measured in 2001, but this was not used for estimating the velocity because it was measured before the first acquisition of the InSAR data. Velocities for site EAY8 were estimated from measurements in 2007, 2009, 2011, and 2015. Before modeling, we reduced the number of InSAR data points by subsampling the two data sets on the basis of the distance from the observed uplift. For the ALOS data, points located within ~10 km of the peak uplift were resampled onto a regular 1 × 1–km grid. Similarly, for the Envisat data, we used all scatterers located within ~10 km of the peak uplift (Fig. 3). We solved for the position, depth, and volume change of the source by performing a hybrid Monte-Carlo downhill simplex inversion, which minimizes the square misfit between the observed and calculated displacements at each scatterer (20). We assumed a Poisson’s ratio of 0.25 and a shear modulus of 30 GPa. Because of the larger number of InSAR observation points, we down-weighed the InSAR observations by a factor of 10 compared to the GPS. We found that additional down-weighing of the InSAR relative (factors of 100 and 1000) to the GPS had little effect on the misfit to the GPS. For the sill, we solved for additional parameters including the position, strike, depth, length, width, and volume change. Given the distribution of data, the extent of the body offshore was difficult to assess. To prevent the model from having unreasonable dimensions, we applied additional penalties if the area of the sill became larger than 500 km2 or had an aspect ratio of more than 3.

Earthquakes in the Matata swarm were relocated using the double-difference hypoDD method of Waldhauser and Ellsworth (41). This approach minimizes the residuals between observed and calculated arrival time differences for pairs of closely spaced earthquakes, using the differential times of catalog P- and S-phase times, as well as differential times derived from waveform cross-correlation of P and S waveforms. We calculated the catalog-based differential times for all event pairs initially separated by less than 8 km and for all seismometers less than 150 km from the cluster centroid using the manually picked P- and S- phase times. Cross-correlation and bispectrum methods were then used to calculate the waveform-based differential times for all event-station pairs after prefiltering, following the technique of Du et al. (42). These derived differential times were weighed on the basis of the quality of the arrival time measurements.

We estimated the spatially varying b value on a ~250-m grid covering the earthquake swarm to have a magnitude of completeness of Mw = 2.3 (fig. S5) using ZMap (43). For each node, a minimum of 50 events above the magnitude of completeness was required to compute a value. We varied the size of the radii around each node but found that for all cases, the largest b values (>1.5) were found closer to the coastline, near our inferred source, and decreased offshore (Fig. 4 and fig. S5).


Supplementary material for this article is available at

fig. S1. Baseline time plot showing the distribution of ALOS and Envisat SAR images (red circles circles) and interferograms (black lines) used in the analysis.

fig. S2. Mean LOS velocities derived using 32 Envisat SAR images of the study region and two offshore islands.

fig. S3. East, north, and vertical displacements recorded at RGMT and RGOP.

fig. S4. Mean LOS velocity derived from ALOS SAR data.

fig. S5. Distribution of model parameters for the Mogi and sill sources.

fig. S6. Cross section through the Matata earthquake swarm along the line A-A′ indicated in the panel on the right of the figure.

fig. S7. Relocated earthquakes from the 2005–2009 Matata earthquake swarm, which are color-coded by year, and available MT solutions for the larger events.

fig. S8. Cumulative number of earthquakes and b value estimates for the Matata earthquake swarm.

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We would like to thank the Japan Aerospace Exploration Agency for access to the ALOS data under project RA4-1093, the ESA for access to the Envisat data under project C1P-14029, and GeoNet ( for the CGPS data. We also thank G. Caldwell and B. Scott for early reviews of the manuscript. We would also like to thank M. Pritchard and an anonymous reviewer for their comments and suggestions, which have greatly improved the manuscript. The aerial photo shown in fig. S2 was sourced from the Land Information New Zealand (LINZ) Data Service ( and licensed by BOPLASS Limited for reuse under the Creative Commons Attribution 3.0 New Zealand license. Funding: This work was supported by public research funding from the Government of New Zealand with additional support from the New Zealand Natural Hazards Research Platform under grant no. 2015-GNS-02-NHRP and LINZ. Author contributions: I.J.H. acquired and processed the InSAR data. S.H. estimated the CGPS velocities. S.B. provided the relocated earthquake catalog. N.P. organized and led the collection of campaign GPS data. I.J.H. carried out the elastic modeling and led the writing of the paper, with contributions from all authors. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.
View Abstract

Stay Connected to Science Advances

Navigate This Article