At the northern Hikurangi margin (North Island, New Zealand), shallow slow slip events (SSEs) frequently accommodate subduction-interface plate motion from landward of the trench to <20 km depth. SSEs may be spatially related to geometrical interface heterogeneity, though kilometer-scale plate-interface roughness imaged by active-source seismic methods is only constrained offshore at <12 km depth. Onshore constraints are comparatively lacking, but we mapped the Hikurangi margin plate interface using receiver functions from data collected by a dense 22 × 10 km array of 49 broadband seismometers. The plate interface manifests as a positive-amplitude conversion (velocity increase with depth) dipping west from 10 to 17 km depth. This interface corroborates relocated earthquake hypocenters, seismic velocity models, and downdip extrapolation of depth-converted two-dimensional active-source lines. Our mapped plate interface has kilometer-amplitude roughness we interpret as oceanic volcanics or seamounts, and is 1–4 km shallower than the regional-scale plate-interface model used in geodetic inversions. Slip during SSEs may thus have different magnitudes and/or distributions than previously thought. We show interface roughness also leads to shear-strength variability, where slip may nucleate in locally weak areas and propagate across areas of low shear-strength gradient. Heterogeneous shear strength throughout the depth range of the northern Hikurangi margin may govern the nature of plate deformation, including the localization of both slow slip and hazardous earthquakes.

Constraints on plate-boundary fault geometry inform slip inversions and studies of seismic hazard and geodynamics. In active subduction zones, the plate interface is usually imaged at a resolution of tens to hundreds of meters using multichannel active-source seismic (MCS) surveys to ~10 km depth (Barker et al., 2009) and at lower (kilometer-scale) resolution and greater depth onshore using passive seismic methods (e.g., Audet and Schaeffer, 2018). We present a detailed onshore three-dimensional (3-D) plate-interface geometry using common conversion point (CCP)–stacked receiver functions (RFs) from a 22 × 10 km array of 49 Guralp 6TD seismometers deployed during the NZ3D experiment (Bell et al., 2018) at the northern Hikurangi margin, New Zealand (Fig. 1). We provide the first detailed 3-D seismic images of a subduction plate interface at 10–17 km depth, deeper than traditional MCS imaging and at higher resolution than most passive methods.

The northern Hikurangi margin represents a significant earthquake and tsunami hazard for North Island, New Zealand (e.g., Shaw et al., 2022). It hosts recurring slow slip events (SSEs) <100 km from the trench (Wallace, 2020), likely associated with fluid migration and subducted seamounts (Shaddox and Schwartz, 2019). Incoming seafloor topography and centimeter-to kilometer-scale protolith variability likely cause locally variable fluid pressures (Chesley et al., 2021), transitional frictional stability (Rabinowitz et al., 2018), geometrical complexity and upper-plate damage (Shaddox and Schwartz, 2019; Sun et al., 2020), and heterogeneous stresses within the deforming zone (Beall et al., 2019), possibly aiding accommodation of the plate-boundary strain budget through aseismic creep and slow slip (Barnes et al., 2020). Direct imaging of the plate interface near the cumulative SSE slip peak allows reevaluation of geodetically inverted slip distributions, fluid content and role in slip, and controls on SSE generation and elastic strain accumulation.

At the northern Hikurangi margin, the Pacific plate subducts beneath the Australian plate at 4–5 cm/yr (Fig. 1; Wallace et al., 2009). Compared to the southern Hikurangi margin, the northern Hikurangi margin currently has low interseismic coupling, a history of Mw ≤7.1 tsunamigenic earthquakes, more frequent shallow M 5 + earthquakes in the past 20 years, and regular SSEs at <20 km depth (Wallace et al., 2009; Wallace, 2020). The ~12-km-thick subducting Hikurangi Plateau (Mochizuki et al., 2021) is geometrically rough; seamounts and volcanic cones comprise porous volcaniclastics and seafloor volcanics atop the oceanic basement and protrude locally through variably thick basin sediments (Barnes et al., 2020). Offshore, low-electrical-resistivity volumes within seamounts locally extend ≤5 km into the Hikurangi Plateau (Chesley et al., 2021). Landward, an underthrust seamount imaged in MCS surveys at 8–10 km depth (Barker et al., 2018) also shows low resistivity (Chesley et al., 2021). The northern Hikurangi plate interface has been imaged to <12 km depth, has locally variable dip (Barker et al., 2009), and locally overlies regions with a high seismic reflectivity, interpreted as fluid-rich sediments (Bell et al., 2010). A regional plate-interface model, extrapolated from offshore seismic data and relocated earthquake hypocenters, dips northwest with dip variations over distances of ≥10 km (Williams et al., 2013). The NZ3D seismometer array was arranged where the plate interface is geodetically unlocked, downdip of the maximum cumulative SSE slip area (Fig. 1) and updip of a localized area of increased coupling (Heise et al., 2017; Wallace, 2020).

We calculated teleseismic RFs using a transdimensional Bayesian framework and a reversible-jump Markov chain Monte Carlo method (Akuhara et al., 2019) on time series bandpass filtered with corner frequencies 0.04 and 5 Hz (≥0.7 km wavelength at 3.5 km/s) and windowed 3 s before and 60 s after the predicted P wave. RFs were depth converted for Ps, PpPs, and PpSs phases by 3-D linear interpolation and ray-path calculation using velocities from Yarce et al. (2021). CCP stacking projected RF amplitudes into a trench-perpendicular volume of 100 × 24 × 90 km with 0.5 × 1 × 0.5 km spacing (Fig. 1). Following Audet and Schaeffer (2018), we averaged and summed amplitudes in each grid cell using phase-weighted stacking, masking contrasting amplitudes between phases; PpPs and PpSs amplitudes were weighted at 3× Ps amplitudes. CCP amplitudes were smoothed perpendicular to the trench using a Gaussian kernel with standard deviation of 2.5 km. We delineated seismic units on 24 trench-perpendicular sections following contrasting or consistent amplitude horizons (Fig. 2A) to form continuous surfaces (Figs. 3A and 3B). Full processing details are provided in the Supplemental Material1.

Surfaces bounding the top and bottom of plate-interface conversions were depth averaged to a mean plate interface (Figs. 2C, 3, and 4G). We evaluated roughness by subtracting mean plate-interface depths from those on a 3-D plane fitted to mean plate-interface depths using a least-squares approach (Fig. 3D). Shear strength of sliding initiation was calculated assuming fault strength is governed by the Mohr-Coulomb failure criterion using isotropic stress at interface depths (Bletery et al. 2017) calculated from estimated density structure with a single friction coefficient of 0.35 (Rabinowitz et al., 2018). Because we concentrate on stress distribution rather than magnitude, we use several convergence azimuths and plunges antithetic to the plate interface, an overburden density of 2300 kg/m3 (Wallace et al., 2019), and pore-fluid factors (fluid pressure/vertical stress) of 0.4, 0.7, and 0.95 (Fig. 4D). Offshore, we used apparent dips on MCS lines (Barker et al., 2009; Gase et al., 2021) and interpolated a margin-perpendicular line from the Williams et al. (2013) model. Strike of the plate interface was assumed to be perpendicular to MCS lines and margin parallel (30°) in the CCP volume. Dip in the CCP volume was calculated on a 5 km rolling average of depth in the margin-normal dimension (parallel to y in Fig. 1) to smooth interface steps produced by CCP resolution.

We interpret the plate interface from a west-dipping positive-amplitude volume (increasing velocity with depth; light red in Fig. 2C) bounded by negative-amplitude volumes. The top (9.5–15 km) and bottom (12–17 km depth) of the plate-interface volume deepen to the northwest from shallow bulges in the east of the study area (Figs. 3A and 3B). Though spacing of the bounding surfaces implies a finite plate-interface thickness, quantifying this and identifying structure within the plate-interface volume requires higher-frequency analysis than used here (conversion amplitudes occur over a 0.5–4 km depth range; see the Supplemental Material text, and Figs. S6–S8), so we use mean depth at each CCP coordinate to describe plate-interface geometry.

Our mean plate interface is 1–4 km shallower than the interface of the Williams et al. (2013) plate-interface model, which does not coincide with any clear horizon within 5–10 km thickness of patchy negative amplitudes in the lower plate (Fig. 2C). This is likely because the Williams et al. (2013) interface was fitted to Vp-sensitive MCS data and seismicity likely concentrated in the lower plate (Reyners et al., 2011), not Vp- and Vs-sensitive RFs which image discontinuities. High Vp/Vs ratios or lower Vs in the wedge would decrease interface depth in our CCP volume. Our mean plate interface bounds the top of relocated seismicity (Reyners et al., 2011; Yarce et al., 2019, 2021; Warren-Smith et al., 2019) and agrees well with the downdip extent of plate-interface depths from recent two-dimensional MCS data (Fig. 2C; Gase et al., 2021). The best-fit plane has a strike/dip of 182°/22°W, ~40° counterclockwise from the strike of the Williams et al. (2013) interface, and is 1–2 km deeper than the mean interface in the east and west and ≤1 km above the mean interface in the study-area center (Fig. 3D). Roughness on the mean interface has wavelengths of 10–20 km and amplitudes of 2–4 km (Fig. 3D); both strike and roughness are consistent with subducting volcanic topography such as subducted seamounts imaged offshore (Bell et al., 2010; Barnes et al., 2020).

The plate interface is thought to exist between siliciclastics and calcareous pelagic sediments in the hanging wall and Hikurangi Plateau volcaniclastics in the footwall (e.g., Barnes et al., 2020), consistent with a velocity increase with depth (Fig. 2). Negative amplitudes below the interface are consistent with low electrical resistivity (Chesley et al., 2021), Vp/Vs >1.83 (Yarce et al., 2021), and low seismic-attenuation (Q) ratios (Qs/Qp <0.75; Nakai et al., 2021) offshore, which suggest the upper 3–5 km of the Hikurangi Plateau hosts fluid-rich locally fractured or porous volcanic basement beneath the plate interface (lowermost boundary in Figs. 2 and 4G). Similar low Qs/Qp ratios in regional data extend downdip across the study area (Fig. 4G), and Vp/Vs >1.83 70–90 km from the trench suggests abundant fluids exist (Eberhart-Phillips et al., 2017; Yarce et al., 2021; Mochizuki et al., 2021).

Wallace et al. (2016) inverted slip distributions at the Hikurangi margin assuming a deeper, smoother plate interface geometry from the Williams et al. (2013) model. The shallower, rougher, mean plate interface presented here may require lower accumulated elastic strain and SSE slip magnitudes in more complex distributions to fit observed geodetic displacements (Wallace et al., 2016; Williams and Wallace, 2018; Perez-Silva et al., 2022), commensurate with greater complexity and slip in fast earthquakes (Shaw et al., 2022). Complex slip distributions would also arise from heterogenous shear strength on a rough plate interface with varied effective normal stress and dip (Fig. 4B; Sun et al., 2020). Interface roughness therefore remains rough to ~16 km landward of the trench, possibly smoothing downdip of the study area with increased shear strain. Frictional parameters may vary more and locally reduce shear strength where roughness controls the lithology hosting sliding, such as in smectite-rich altered basalts in the lower plate (Wallace et al., 2019; Barnes et al., 2020). Slip on faults with heterogeneous shear strength and frictional slip behaviors nucleates in, and propagates from, the weakest or highest-stress areas (Fig. 4; Luo and Ampuero, 2018; Beall et al., 2019). Our variable geometry could facilitate slow slip by limiting slip nucleation and acceleration locally to weak or high-stress areas. Heterogeneous strength and (locally) elevated pore-fluid factors (e.g., Luo and Ampuero, 2018) mean multiple distinct areas are likely close to failure (Figs. 4B and 4D), allowing sequential local slip nucleation in several areas or loading from slip elsewhere to trigger slip across large areas of the plate interface (Perez-Silva et al., 2022).

The 2014 CE Gisborne SSE initiated atop a subducted seamount 20–35 km from the trench and propagated downdip, where the greatest slip magnitudes occurred (Fig. 4; Wallace et al., 2016; Yohler et al., 2019). Maximum slip magnitudes occurred where low megathrust dip causes low shear-strength gradients normal to the trench, conducive to slip propagation (see the Supplemental Material text; Fig. 4; Bletery et al., 2017). A shallower plate interface would have less overburden and lower effective normal stress than the Williams et al. (2013) model, making seismic and/or slow slip weaker and more sensitive to pore-fluid pressure changes reducing the range in shear stress (Fig. 4D; Shaddox and Schwartz, 2019; Warren-Smith et al., 2019). We speculate that variable shear strength from plate roughness may link to variable fault slip style along this margin. High-resolution studies are needed across plate boundaries that do and do not experience slow slip to test correlations between roughness and slip style.

Instruments were loaned by SEIS-UK (Leicester, UK; loan 1039), funded by UK Natural Environment Research Council (NERC) grant NE/MO21203/1. H. Leah and Å. Fagereng are supported by European Research Council Starting Grant 715836 (“MICA”); R. Bell is supported by NERC grant NE/MO21203/1. B. Fry, S. Henrys, and K. Jacobs are supported by the New Zealand Ministry of Business, Innovation and Employment Strategic Science Investment Funds and Endeavor Fund grant (Hikurangi subduction earthquakes and slip behavior and C05X1605) to GNS Science. We thank all those involved in the NZ3D-FWI (New Zealand 3-D full waveform inversion) data collection, including M. Cave, the Gisborne District Council, and the Gisborne landowners who hosted seismic stations. Seismograms were assessed using RF and ObsPy (Beyreuther et al., 2010; Eulenfeld, 2020; Krischer et al., 2015). We thank P. Barnes and two anonymous reviewers for their comments, which helped improve the manuscript.

1Supplemental Material. Details of seismic arrivals, RF processing, CCP stacking, synthetic CCP tests, and shear strength modeling. Please visit https://doi.org/10.1130/GEOL.S.20465385 to access the supplemental material, and contact [email protected] with any questions.