Seismic Constraint on Heterogeneous Deformation and Stress State in the Forearc of the Hikurangi Subduction Zone, New Zealand Open Access

The Hikurangi subduction zone (HSZ) spans nearly 500 km along the eastern coast of the North Island of New Zealand, and marks the collisional boundary between the Pacific and Australian tectonic plates (Fig. 1). The Pacific plate subducts transpressively westward beneath the North Island, at the southern end of the Tonga–Kermadec–Hikurangi subduction system. The section of the Pacific plate subducting along the HSZ consists of a Cretaceous large igneous province (the Hikurangi plateau) with a progressively southward thickening (10–15 km thick) buoyant oceanic crust (Mochizuki et al., 2019). This ultimately resists subduction and causes a complex southward transition into the strike‐slip Alpine fault across the South Island of New Zealand.

The interface (megathrust fault) between the Pacific and Australian plates at the HSZ is believed to be capable of producing magnitude 9.0 earthquakes and associated tsunamis, with geological evidence of four megathrust ruptures over the past 2000 yr (Pizer et al., 2021). However, Global Navigation Satellite Systems data over the North Island suggest that only the southern part of the megathrust is currently locked and accumulating elastic strain that can be released in a large earthquake (Wallace et al., 2004). In contrast, the Australian and Pacific plates are believed to be uncoupled, and the megathrust creeping, in the northern part of the HSZ (Fig. 1). Furthermore, slow‐slip events (SSEs) on the Hikurangi megathrust fault have been observed (Wallace, 2020) following the pattern of plate locking along the margin but generally down‐dip of the locked part of the fault (Fig. 1). These observations are consistent with modeling that suggests SSEs occur in a transitional frictional state (Liu, 2013). Changes in the location of SSEs as well as the inferred state of plate coupling along the margin indicate that the mechanical conditions within the subduction zone change significantly from north to south (Wallace, 2020).

Knowledge of the structural and mechanical conditions of currently locked and creeping segments of megathrust faults is important for understanding the seismic hazards posed by these tectonic systems. This research uses passive seismic data to constrain the tectonic grain over the forearc of the HSZ to address how present‐day and inherited deformation are distributed within the context of the spectrum of megathrust behavior along the margin. Specifically, this work examines the directional dependence of seismic scattering structure, and explores the spatial variability of these data in relation to forearc deformation, plate coupling, and subduction geometry and morphology.

Distant earthquakes produce teleseismic P waves that travel through the Earth’s mantle and arrive almost vertically beneath a seismograph station. Some of the P‐wave energy will be scattered by structure beneath the seismograph station and will be converted into S‐wave phases that appear on the horizontal‐component channels. Deconvolving the vertical component from the two horizontal components of motion isolates these converted S waves (by approximating the vertical component as a source wavelet). Hence, a three‐component teleseismic P waveform contains a record of the structure beneath a seismic station in the form of a P–S receiver function (RF) (e.g., Langston, 1979). A single RF is a waveform that approximates the impulse response of the structure beneath a seismic station. The timing of a particular pulse in the RF waveform constrains the depth to the structure that generated the pulse. The amplitude of the pulse constrains the magnitude of the velocity contrast across the structure.

In this work, we compute RFs for 53 long‐running seismic stations (Fig. 1) part of the national seismograph network of New Zealand (GeoNet) using recordings of earthquakes with magnitudes > 5.5 that occurred between January 2000 and March 2020, within epicentral distances of 30°–90° for P arrivals and 90°–150° for PP arrivals, relative to each station. For each earthquake, we use the IASP91 Earth model (Kennett and Engdahl, 1991) to extract 120 s time windows centered on the predicted arrival (P or PP). The use of PP arrivals is not conventional in RF analysis, because these phases typically have lower amplitudes. We consider these phases to improve azimuthal data coverage at each station. However, due to the quality control measures used in this study (discussed subsequently), there are 2–3 orders of magnitude more P data than PP data in our final dataset. Examples of P and PP RFs are shown in the supplemental material, available to this article. Many of the GeoNet stations used in this work are short‐period instruments, which are typically inferior to broadband instruments in RF applications. However, many recent studies have successfully applied RF analysis using similar instruments (e.g., Niu et al., 2005).

All data time windows are initially corrected for their instrument response. For each event, the signal‐to‐noise ratio (SNR) is calculated as the ratio of root mean squared amplitudes on the vertical channel, and compares the 30 s of data before and after the predicted P (or PP) arrival (measured over the frequency band of 0.05–1 Hz). Given the long recording times available for the instruments considered in this work, we only consider events with SNR greater than 5 dB. Horizontal components are then rotated into radial and transverse orientations, relative to the event location. To isolate P–S converted phases (i.e., compute RFs), we use a modified Wiener spectral deconvolution that is regularized by the pre‐event covariant noise spectrum between vertical and horizontal components. This has the advantage of minimizing potential contamination of the RFs due to seasonal noise effects (Audet, 2010). RFs are then band‐pass filtered using a second‐order Butterworth filter between 0.05 and 0.75 Hz. Examples of radial and transverse RFs computed at station MRHZ are shown in Figure 2, sorted by event back azimuth (BAZ).

For isotropic media with horizontal velocity discontinuities, all of the teleseismic P‐waveform energy is contained within the radial plane (i.e., no energy appears on transverse RFs). However, in the presence of material anisotropy or dipping structure, some of the teleseismic P‐waveform energy can be refracted out of the radial plane and onto the transverse plane. Consequently, the observation of coherent pulses on transverse‐component RFs is indicative of anisotropy and/or dipping structure beneath a seismic station (e.g., Savage, 1998). For all stations considered in this work, we observe coherent pulses on transverse‐component RFs at P–S delay times (i.e., depths) associated with the crust of the overriding Australian plate (Fig. 2).

In the presence of anisotropy or dipping structure, the amplitude and timing of pulses on both the radial and transverse RFs depends on the orientation of the structure (i.e., strike of dipping structure or fast‐axis direction of anisotropy) relative to the BAZ of the event used to compute the RF. Forward modeling of RFs suggests that dipping structure or a plunging axis of anisotropy characterized by hexagonal symmetry produces pulses that vary in amplitude and timing with a 360° periodicity with respect to BAZ (e.g., Savage, 1998; Porter et al., 2011; Schulte‐Pelkum and Mahan, 2014). Such pulses are observed in our RF dataset (Fig. 2). A heuristic approach to constrain the orientation of the structure beneath a seismic station, without requiring complex waveform modeling, is to decompose the RF data at a given station into a set of harmonics with respect to BAZ (e.g., Bianchi et al., 2010; Schulte‐Pelkum and Mahan, 2014; Audet, 2015; Schulte‐Pelkum et al., 2020). Previous studies considering RF data for stations throughout the HSZ have identified significant and complex seismic scattering structure in the overriding plate, further motivating the use of such an approach (e.g., Savage, 1998; Bannister et al., 2007; Savage et al., 2007).

In a geographic reference frame, the coordinate system of the BAZ harmonics is defined by sine and cosine functions that are aligned with the north–south and east–west directions. This implies that dipping structure (or plunging anisotropy) that is aligned in a north–south orientation will be described by only one of the B(t) components (i.e., the complementary B(t) component contains no signal). Structure that is not aligned in either a north–south or east–west direction will produce a signal on both B(t) components. By introducing the phase shift $α$ in equation (1), the reference frame can be rotated to determine the orientation of the coordinate system that produces maximum signal amplitude on a particular B(t) component (Audet, 2015) and, therefore, the orientation of geologic structure. In this study, we search for the angle $α$ that maximizes signal variance on the $B2(t)$ component to estimate the strike orientations of dipping structure or the trends of the symmetry axes of plunging anisotropy. We specifically analyze RF data over the time window from 0.5 s (to avoid the direct P arrival) to the P–S delay time associated with the depth of the plate interface based on the model by Williams et al. (2013), which is unique beneath each station. The P–S delay times are estimated using a unique velocity profile for each station, based on a 3D seismic tomography model of New Zealand (Eberhart‐Phillips et al., 2020). Consequently, our results represent the dominant orientation of structure within the crust of the overriding Australian plate. A map of predicted P–S delay times for the region is shown in the supplemental material.

We calculate harmonic decompositions of RF data for all stations over the forearc of the HSZ. The A(t) harmonic components for all stations, which represent the isotropic structure beneath the stations, are shown in the supplemental material. For each station, we estimate the orientation of the strike of dipping structure, or trend of plunging anisotropy, within the overriding Australian plate (Fig. 3a). To determine the uncertainty in the estimated orientations, we perform a bootstrap analysis by randomly resampling the RF data for each station (with replacement) and recomputing the harmonic decompositions 100 times. The black bow ties in Figure 3a represent the 80% confidence interval on the estimated orientation for each station. Overall, our results show good agreement between proximal stations throughout the study region. Throughout the southern HSZ, our estimated orientations are generally aligned with the margin as well as the strike of the plate interface shown as contour lines in Figure 1 based on the model by Williams et al. (2013). However, the orientations of structure beneath stations in the northern HSZ are predominantly aligned in a west‐southwest–east‐northeast direction.

Some studies discriminate the signature of a dipping interface (as opposed to plunging anisotropy) based on the presence of an azimuthally periodic phase on the transverse RFs at zero delay time (e.g., Porter et al., 2011). This work investigates structure in a subduction zone, where dipping structures are expected below the overriding plate. Using an initial time‐window bound of 0.5 s, we avoid analyzing this (possible) phase, which is likely influenced by deeper structure that is not considered here. It is also possible (in some cases) to discriminate dipping structure from plunging anisotropy based on associated phases on the A(t) component, as anisotropy may not produce such energy (Licciardi and Piana Agostinetti, 2016). However, the heuristic approach used in this work considers the harmonic signal in the RF data over a broad time window, rather than correlating individual phases observed at the same delay time but on different harmonic components. For these reasons, we are not able to distinguish between plunging anisotropy and dipping structure from this analysis alone.

Forward modeling of RF data suggests that the signature for anisotropy with a plunging slow axis of symmetry is similar to the signature of anisotropy with a plunging fast axis of symmetry (with trend and plunge rotated 180° and 90°, respectively). With real data, it is not reasonable to expect to discriminate between these two cases, resulting in an ambiguity of 180°. Furthermore, as this analysis estimates $α$ by maximizing signal variance on the $B2(t)$ component, the equivalent signal variance is found with $α$ rotated by 180°. Finally, the estimated orientations represent the alignment of the dominant structure (or structures) beneath each station within the overriding plate. Hence, our analysis is not able to distinguish multiple dipping structures or multiple anisotropic layers (or combinations thereof).

Figure 3b shows the magnitudes of the $B2(t)$ harmonic components for each station. Specifically, these values are the standard deviations of the $B2(t)$ components over the specified time windows corresponding to depths of the overriding Australian plate (i.e., from 0.5 s to the P–S delay time associated with the plate interface, as discussed earlier). These values are normalized by the total signal from all harmonic components to represent the percentage of the signal in the RF dataset (at a given station) that is described by a 360° periodicity with respect to BAZ. Hence, this describes the relative strength of the signal associated with dipping structure or plunging anisotropy. For stations in the northern HSZ, approximately 25%–30% of the RF data variance can be described by dipping structure or plunging anisotropy. In the northern part of the southern HSZ, over 40% of the RF dataset can be described by dipping structure or plunging anisotropy, decreasing to approximately 20% toward the southernmost point of the North Island.

The estimated orientations, and relative magnitudes, from the harmonic decompositions of RF data suggest significant margin‐parallel heterogeneity in tectonic fabric within the forearc. Previous forward modeling of RF data suggests that the magnitude of the $B2(t)$ harmonic component is influenced by the magnitude of the velocity contrast across a dipping interface, the magnitude of plunging anisotropy, as well as the dip or plunge angle of such structures (e.g., Porter et al., 2011; Schulte‐Pelkum and Mahan, 2014). In our case, the strength of the $B2(t)$ harmonic signal may also be influenced by the superposition of multiple structures. It is unlikely that the variability in our results over the entire forearc can be explained by a single factor. Several recent studies observe along‐margin variability in seismic velocities (e.g., Reyners and Eberhart‐Phillips, 2009; Eberhart‐Phillips et al., 2010), seismic attenuation (e.g., Eberhart‐Phillips et al., 2017), and electrical resistivity (e.g., Heise et al., 2013), suggesting our estimated changes in the relative magnitude of the $B2(t)$ harmonic components may be partially attributed to heterogeneity in bulk geophysical properties along the margin. However, based on our estimated structural orientations, we suggest that the dominant source of the $B2(t)$ harmonic signal is different between the southern HSZ and the northern HSZ, and that observed changes across the southern HSZ (where the plate interface is coupled) likely reflect changes in subduction morphology.

Motion of the Pacific plate relative to the Australian plate is oblique to the margin (Fig. 1). This obliquity increases southward due to the changing orientation of the margin as well as a progressive rotation of relative plate motion (Kreemer et al., 2003). Margin‐parallel relative plate motion is accommodated predominantly by a general clockwise rotation of the entire forearc (Wallace et al., 2004), as well as by a belt of major dextral strike‐slip faults spanning the entire forearc (Litchfield et al., 2014), as shown in Figure 3. Motion on these faults is believed to accommodate a significant portion of margin‐parallel relative plate motion in the southern HSZ, decreasing to almost no motion in the northern HSZ (Litchfield et al., 2014). Active‐source seismic imaging suggests that these faults are splayed and likely intersect the plate interface at depth (e.g., Henrys et al., 2013). Hence, these faults are expected to be dipping structures, and shearing along these faults could impart a foliation within the crust of the overriding Australian plate, generally aligned with their strike orientation (which is similar to the strike of the plate interface in this region).

We suggest that the observed margin‐parallel variations in $B2(t)$ harmonic orientations and magnitudes in the southern HSZ, where the plate interface is coupled, are due to changes in the morphology of the plate interface as well as a progressive change toward transpressional deformation southward (Fig. 4). This is consistent with results from Graham et al. (2020) who map anisotropy fast‐axis orientations over the southern HSZ using measurements of shear‐wave splitting from local earthquakes (via splitting tomography). Their results generally show margin‐parallel orientations, with a deviation at the southern tip of the North Island that suggests a change in the source of anisotropy. We also note that the southward decrease in $B2(t)$ harmonic magnitudes in the southern HSZ (Figs. 3b and 4b) correlates with a significant change in alignment of the deformation front offshore and corresponding changes in plate interface strike (Fig. 1). Consequently, the geometry and morphology of deformation within the overriding plate is expected to vary accordingly across this region, altering the geometry and potentially the strength of crustal foliations. However, long time periods may be required to impart foliations within the forearc crust, and the orientations of structures relative to the strike of the plate interface (and offshore deformation front) are not constant throughout the geologic history of the HSZ.

In the northern HSZ, where the plate interface is not coupled, we suggest that the orientations estimated from the harmonic decomposition of RF data are not related to subduction directly. Instead, we propose that the general clockwise rotation of orientations toward an west‐southwest–east‐northeast direction may be related to extension and rifting of the Taupo Volcanic Zone (TVZ, Fig. 3a) further inland. We note that our results are consistent with spatially averaged estimates of fast‐axis orientations of anisotropy derived from measurements of shear‐wave splitting from local earthquakes (via splitting tomography), which reveal a rotation in structural orientations along the margin toward a similar west‐southwest–east‐northeast direction in the northern HSZ (Illsley‐Kemp et al., 2019). The orientations may represent structural fabrics imparted by past deformation that have subsequently been rotated clockwise with the overall rotation of the forearc (Wallace et al., 2004). The rotation of orientations toward an east–west alignment for the northernmost stations on the Raukumara Peninsula (northernmost stations in Fig. 3a) may be related to changes in tectonic stability of the peninsula relative to the rest of the Australian plate, as supported by observed changes in paleomagnetic data (Mumme et al., 1989) and earthquake focal mechanisms (Reyners and McGinty, 1999) over this region.

The observed orientations in the northern HSZ may alternatively be related to present deformation, rather than structural fabrics imparted by past deformation. Extension in a north‐northwest–south‐southeast direction suggests that the direction of maximum horizontal stress is oriented in an orthogonal (i.e., west‐southwest–east‐northeast) direction, consistent with our observations. In this case, our results may represent the preferential opening of fluid‐filled microfractures that are oriented in a west‐southwest–east‐northeast direction within the extensional forearc crust. This would result in a symmetry axis of anisotropy consistent with our observations. Furthermore, geochemical analysis of onshore spring waters throughout the forearc suggests greater permeability of the upper plate in the northern HSZ compared to the southern HSZ, which may be related to changes in tectonic stress state (i.e., compression vs. extension) along the margin (Reyes et al., 2010; Barnes et al., 2019). Finally, changes in stress state offer a possible explanation for the aforementioned observations of along‐strike variability in bulk geophysical properties (e.g., Reyners and Eberhart‐Phillips, 2009; Heise et al., 2013; Eberhart‐Phillips et al., 2017), which may reflect changes in signature of crustal fluids.

The HSZ displays enigmatic margin‐parallel changes in plate coupling and SSE locations, suggesting heterogeneity in the physical and frictional state of the megathrust. In this work, we compute RFs from teleseismic data recorded at 53 long‐running seismic stations over the forearc of the HSZ to constrain seismic scattering structure. We observe a strong BAZ dependence in the timing, amplitude, and polarity of RF pulses at P–S delay times associated with depths of the overriding forearc crust. We decompose the RF data into harmonic components to constrain the orientations of associated structures that cause these signals. We observe significant along‐strike variability in the magnitudes and orientations of $B2(t)$ harmonic signals, suggesting heterogeneous deformation throughout the forearc. In the southern HSZ, we suggest that our observations are related to splayed dextral faults cutting through the forearc crust, which change in morphology along strike. In the northern HSZ, where the overriding forearc crust is uncoupled from the subducting Pacific slab, we suggest that our observations are related to past deformation, clockwise rotation, or extension of the forearc crust. Consequently, our results provide further evidence of variability in tectonic stress state of the forearc crust along the margin, which is likely related to the state of interseismic coupling. This work yields further evidence of the heterogeneous structure and dynamics of the HSZ, and provides useful constraints on the relationships among megathrust coupling, forearc deformation, and seismic hazards.

GeoNet seismic data are freely available at www.geonet.org.nz and were downloaded using ObsPy packages (Beyreuther et al., 2010). Receiver functions (RFs) were processed using the RFPy package (Audet, 2020). Analysis results (⁠$B2$ magnitudes and orientations) with associated uncertainties are listed in the supplemental material.

The authors acknowledge that there are no conflicts of interest recorded.

This work is supported in part by the Natural Science and Engineering Research Council (Canada) through a Vanier Canada Graduate Scholarship to Jeremy M. Gosselin and a Discovery Grant to Pascal Audet. This work is also supported by the R‐CET Endeavour program, the Hikurangi Subduction Margin Endeavour Fund, and the GNS Strategic Science Investment Fund. The authors thank Martha Savage and an anonymous reviewer for useful comments that improved this article.