The northern Hikurangi margin three-dimensional plate interface in New Zealand remains rough 100 km from the trench

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 earth-quake 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.


INTRODUCTION
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 centimeterto 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 aseis-mic 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.

TECTONIC SETTING
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 M w ≤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 fluidrich 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).

METHODS
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 trenchperpendicular 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-perpen-dicular 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/m 3 (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.

PLATE INTERFACE GEOMETRY
We interpret the plate interface from a westdipping 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 plateinterface 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., 2019Yarce et al., , 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).

MECHANICAL IMPLICATIONS OF A SHALLOWER, ROUGH INTERFACE
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.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.
Figure 3. (A-C) Contoured plate-interface geometry.(D) Depth from the mean plate interface to the best-fit plane.The mean plate interface is shallower in the red areas and deeper in the blue areas.Fitted plane strike and dip (182°/22°W) is shown in the lower left.Inverted triangles are NZ3D seismic-network stations.

Figure 4 .
Figure 4. Plate interface depth and mechanics across the Hikurangi margin (New Zealand).(A) Plate interface depth from this study and SHIRE (Gase et al., 2021) and 05CM-04 (Barnes et al., 2020) multichannel active-source seismic lines.Gray fill is the Gisborne 2014 CE slow slip event (SSE) slip initiation area (Yohler et al., 2019); red fill is the subducted seamount (Bell et al., 2010).(B) Shear strength (τ) of initial sliding at hydrostatic pore-fluid pressure on the mean plate interface (see text for details).(C) SSE slip magnitudes for the 2014 Gisborne event and for 2002-2014 (Wallace et al., 2016; Wallace, 2020) along transect X-X′ shown in A. (D,E) Shear strength of frictional sliding along transect shown in A for imaged plate interface (solid lines near trench are from 05CM-04 line; dashed lines are from SHIRE line; solid lines onshore are from this study) and Williams et al. (2013) model (dotted line) at various pore-fluid factors (λ) (D) and various stress orientations (E).Colored fill in D shows margin-parallel range on mean plate interface.(F) Shear-strength gradient for frictional sliding with various stress orientations (key in E). (G) Geometry of plate interface and geophysical data along transect shown in A. Bathymetry (gray line) and plate interfaces of Williams et al. (2013; red dotted line), Gase et al. (2021; red dashed line) within 5 km laterally of the section line, Barnes et al. (2020; solid red line near trench), and this study (solid red line onshore) are plotted alongside base of negative amplitudes in lower plate (black line; Fig. 2) and margin-parallel upper-and lowerbound depth ranges of mean plate interface and base of underlying negative amplitudes (black outline and gray fill, respectively).Green dashed line is the base of the Hikurangi Basement (unit HKB of Barnes et al., 2020).Gray dots are seismicity within 10 km of the section line (Reyners et al., 2011; Yarce et al., 2019, 2021; Warren-Smith et al., 2019).Grey vertical bar 20-35 km from the trench in C-G is the slip initiation area of the Gisborne 2014 SSE (Yohler et al., 2019).