Three‐Dimensional P Wave Velocity Structure of the Northern Hikurangi Margin From the NZ3D Experiment: Evidence for Fault‐Bound Anisotropy

We present a high‐resolution three‐dimensional (3‐D) anisotropic P wave velocity (Vp) model in the northern Hikurangi margin offshore Gisborne, New Zealand, constructed by tomographic inversion of over 430,000 first arrivals recorded by a dense grid of ocean bottom seismometers. Since the study area covers a region where shallow slow slip events (SSEs) occur repeatedly and the subduction of a seamount is proposed, it offers an ideal location to link our understanding of structural and hydrogeologic properties at megathrust faults to slip behavior. The Vp model reveals an ~30‐km‐wide, low‐velocity accretionary wedge at the toe of the overriding plate, where previous seismic reflection studies show a series of active thrust faults branching from the plate interface. We find some locations with significant Vp azimuthal anisotropy >5% near the branching faults and the deformation front. This finding suggests that the anisotropy is not ubiquitous and homogeneous within the overriding plate, but more localized in the vicinity of active thrust faults. The fast axes of Vp within the accretionary wedge are mostly oriented to the plate convergence direction, which is interpreted as preferentially oriented cracks in a compressional stress regime associated with plate subduction. We find that the magnitudes of anisotropy are roughly equivalent to values found at oceanic spreading centers, where the extensional stress regime is dominant and the crack density is expected to be higher than subduction zones. This consideration may indicate that additional effects such as fault foliation and clay mineral alignment also contribute to upper plate anisotropy along subduction margins.


Introduction
The Hikurangi subduction margin, New Zealand (Figure 1), has not experienced any major (Mw > 7.2) subduction interface earthquakes since historical record keeping began~160 years ago (Webb & Anderson, 1998). The overriding Australian plate, by contrast, has produced 10 large (Mw > 7.0), damaging earthquakes over the same period (Webb & Anderson, 1998). Paleoseismic interpretations reveal that up to 10 subduction earthquakes may have occurred in the past 7,000 years, with the strongest evidence for a full Hikurangi margin rupture occurring 870-815 years before present, potentially with Mw > 9.0 (Clark et al., 2019). However, a diverse array of aseismic creep and slow slip event (SSE) behavior has also been documented along the Hikurangi subduction interface (Wallace, 2020). In the northern part of the margin, offshore Gisborne, frequent (18-to 24-months recurrence intervals), small (moment releases of Mw = 6.3-6.8) SSEs occur at <10-to 15-km depths Wallace et al., 2016). This area is also known to have hosted two large tsunami earthquakes in 1947 (Bell et al., 2014;Doser & Webb, 2003). Despite the importance for assessing seismogenic and tsunami potentials, the physical properties of the shallow megathrust faults and splay faults and their three-dimensional variation are poorly documented and underlying mechanisms controlling the spatial variability in megathrust slip behavior remain largely unknown. Figure 1. (a) Regional tectonic map of New Zealand. The black solid line is the plate boundary between the Pacific plate and Australian plate from Bird (2003). The black box indicates the location of 1b and 1c. (b) Seismic activity in the study area. Black dashed contours indicate the slip amount (in millimeter) of the SSE in 2014 (Wallace et al., 2016). Yellow dots, red stars, and a green triangle are the locations of tremor , burst-type repeating earthquakes (Shaddox & Schwartz, 2019) and the tsunami earthquake in 1947 (Doser & Webb, 2003), respectively. The location of inferred subducting seamount (orange shade delineated by the red line) is from Barker et al. (2018). The black arrow shows the convergence direction of the Pacific Plate relative to the North Island at the Hikurangi Trough (Wallace et al., 2004). (c) Layout of the NZ3D OBS experiment. White circles indicate the locations of OBSs. Numbered locations correspond with OBSs indicated in Figure 2. White circles with a colored circumference (red, purple, and blue) correspond to three groups on the difference geological units (incoming plate, accretionary wedge, and forearc basin) that are used to show the spatial variation in traveltime residual in Figure 6. Dark green lines are locations of air gun sources. Yellow squares are the locations of four IODP Exp372/375 drill sites . Brown line shows the multichannel seismic reflection profile of 05CM-04 (Bell et al., 2010;Barker et al., 2018;Gray et al., 2019;Figure 1d). The red dashed line marks the outline of location of the inferred subducting seamount . (d) The 05CM-04 seismic reflection image Wallace et al., 2019).
In 2017-2018, an intensive seismic imaging project NZ3D was carried out in the northern Hikurangi margin (Figure 1c). This project consisted of three components; a 3-D multichannel seismic (MCS) reflection profiling (MGL1801 participants, 2018), a 3-D high-resolution ocean bottom seismograph (OBS) experiment (Kellett et al., 2018), and a 3-D onshore seismic observation . The NZ3D aims to comprehensively understand the 3-D structural, stratigraphic, and hydrogeologic properties along the subduction megathrust faults where a range of slip behavior is observed. The offshore NZ3D survey area overlaps the previously acquired 2-D MCS reflection profile of 05CM-04 Bell et al., 2010; Figure 1), ongoing passive seismic and geodetic monitoring studies (Shaddox & Schwartz, 2019;Todd et al., 2018;Wallace et al., 2016) and International Ocean Discovery Program (IODP) Drilling Expeditions 372 and 375 Wallace et al., 2019). The offshore Gisborne area thus provides one of the most comprehensive geophysical/geological datasets across a subduction zone and enables understanding of properties and behaviors of shallow subduction faults. In this paper, we present a detailed P wave velocity (Vp) model derived from the NZ3D OBS experiment. We determine compressional wave azimuthal anisotropy, which contributes to understanding crustal stress heterogeneity at lateral scales of a few kilometers. Our results are compared with logging data, borehole breakouts, shear wave splitting anisotropy, and mapped faults. This comparison is important since this may constrain possible variations in structure, stress, and their interpretation with depth at subduction zones.

Geological Setting
Along the~500-km-long Hikurangi subduction zone, the oceanic Hikurangi Plateau, a large igneous province of Cretaceous age, is subducting beneath the North Island of New Zealand (e.g., Wood & Davy, 1994). The Hikurangi Plateau includes numerous seamounts and other volcanic features and is composed of basaltic basement thickly overlain by sedimentary layers of Mesozoic to Cenozoic ages (Barnes et al., 2020;Davy et al., 2008). Recent seismic surveys confirmed that the crust is over 10 km thick, and its thickness gradually increases from north to south in accordance with the latitudinal variation in seafloor depth Mochizuki et al., 2019). In the northern part of the plateau where the seafloor is deeper (>3,000 m), more seamounts emerge on the seafloor, but the southern part of the plateau also hosts a number of basement highs that are buried by thick sediments Mochizuki et al., 2019). Correspondingly, the thickness of the trench-fill sediments varies significantly along the trench from approximately 1 km in the north to >5 km in the south, which drives a progressive increase in the magnitude and extent of frontal accretion along the margin (Barnes et al., 2010).
The subduction of the Hikurangi Plateau at the Hikurangi margin is thought to have been initiated in the late Oligocene to early Miocene (Ballance, 1976;Lamb & Bibby, 1989;Rait et al., 1991). The current convergence rate progressively decreases from 58 mm/year in the north to 19 mm/year in the south (Wallace et al., 2004). The direction of the convergence is almost perpendicular to the trench in the north, but it is highly oblique to the trench strike in the south. The present interplate coupling rate also changes latitudinally with a sudden transition at~40°S: To the north, the margin is characterized by low interplate coupling and shallow (<10 km) SSEs. To the south, a zone of strong coupling extends to depths of >40 km and is flanked downdip by regions hosting SSEs (Wallace, 2020). Seismic reflection studies recognized several subducting seamounts and high reflectivity zones in the source regions of the SSEs (Barker et al., 2009;Bell et al., 2010). The reflective zones are typically interpreted to be fluid-rich sediments entrained deep into the subduction zone by the subducting seamounts.
In the northern margin where the plate interface experiences SSEs at shallow depth, the forearc is characterized by a narrow-imbricated structure consisting of late Cenozoic accretionary prism at the toe of the overriding plate and landward Cretaceous to Paleogene rocks, both of which are variably covered by Miocene to recent sedimentary units (Barnes et al., 2010). The study area of our 3-D seismic experiment covers the Tuaheni ridge and the northern Tuaheni landslide complex on the upper continental slope. Within the accretionary prism, active landward-dipping thrust faults that splay off the plate interface are recognized (Bell et al., 2010;Barker et al., 2018;Figure 1d). Multichannel seismic (MCS) reflection and seafloor magnetic studies by Bell et al. (2010) and Barker et al. (2018) suggest that a large (30 × 10 km) seamount is subducting in this region and shallow tectonic tremors detected by Todd et al. (2018) are mostly located just above and landward of the suspected seamount ( Figure 1c). Using the same reflection data as Barker et al. (2018), Gray et al. (2019) derived a fine-scale velocity structure of the uppermost 1.5-2 km below the seafloor by applying a full-waveform inversion (FWI) technique and revealed sudden velocity changes within the frontal accretionary wedge indicating low-velocity zones along the splay faults, which may be related to fluid flow.

Seismic Anisotropy in the Crust of the Northern Hikurangi Margin
Stress field parameters can be measured using borehole breakouts, focal mechanisms of local earthquakes, shear wave splitting, and compressional wave azimuthal anisotropy. In the crust, anisotropy is mainly caused by either preexisting structural features, such as fault zones, or by the preferential closure of fluid filled fractures in response to the direction of horizontal compressive stress (Crampin, 1981). Differential horizontal stresses preferentially open the cracks that are oriented parallel to the maximum horizontal compressive stress (SH max ). Consequently, the magnitude of crack-controlled seismic anisotropy can be directly related to the magnitude of SH max (e.g., Zatsepin & Crampin, 1997). It is also known that seismologically measurable anisotropy can be caused by the orientation of minerals in the rocks even in the absence of cracks. For example, Okaya et al. (1995) used rock samples from the Alpine Fault Zone in New Zealand to show that Haast schists and Torlesse greywackes exhibit several percent of anisotropy in Vp even in the plane of the foliation. Sedimentary rocks that are enriched in clay minerals and form a fine layering are another major source of seismic anisotropy at shallow depths (Wang, 2002).
The northern Hikurangi margin also provides an opportunity to examine the present-day variation in stress state in the upper and downgoing plates, knowledge of which is critical for understanding the storing and release of seismic energy (Townend et al., 2012;Zoback, 1992). Regional shear wave splitting analyses along the northern Hikurangi margin using local earthquakes show mainly NE-SW trench-parallel fast azimuths on the North Island (Audoine et al., 2004;Townend et al., 2012;Unglert, 2011). The NE-SW orientation of SH max is also observed in the subducting plate beneath the North Island by Warren-Smith et al. (2019), using earthquake focal mechanisms, which is consistent with a normal faulting stress regime in the slab. By contrast, a comprehensive catalogue of anisotropy measurements based on passive seismic data along the margin (Illsley-Kemp et al., 2019), sampling the upper plate forearc, observe predominantly plate-motion parallel (E-W) fast orientations in the northern Hikurangi margin. Eberhart-Phillips and Reyners (2009) found that the Hikurangi margin forearc is characterized by heterogeneous anisotropy (both in magnitude and azimuth) using earthquake arrival time data. Shear wave splitting for event-station pairs on HOBITSS ocean bottom seismometers and onshore GeoNet seismic stations yields E-W fast azimuths in the accretionary wedge close to the trench and variable directions further landward, which may be explained by fracture and fault patterns created by subducting seamounts (Zal et al., 2020). Borehole breakouts at drill sites support that SH max at shallow depths in the accretionary wedge is orthogonal to the deformation front .
Shear wave splitting suffers poor spatial resolution in the depth direction since splitting properties are cumulative along the path for delay time. Shear waves may be re-split if they encounter a later and differing anisotropic layer, then a shallow secondary polarization instead of a deeper primary may be measured. Although a 3-D seismic refraction survey typically can sample only upper several kilometers crust, this method can overcome the issue of depth resolution by quantifying the variation in Vp as a function of the azimuth and depth of seismic raypaths throughout the 3-D volume. This approach can provide an accurate estimate of compressional wave azimuthal anisotropy at shallow depths (e.g., Dunn, 2015).

Seismic Data Acquisition and Description
The NZ3D OBS experiment acquired 3-D seismic refraction/wide-angle reflection data with multiazimuth ray coverage using the largest number of OBSs and densest air gun shooting in subduction zones to date. The seismic data were collected via a series of three marine expeditions. In December 2017, 100 OBSs of JAMSTEC were deployed by R/V Tangaroa with an average spacing of 2 km on four parallel lines across the Hikurangi Trough ( Figure 1c). Each OBS contained hydrophone and three-component geophone sensors with natural frequency of 4.5 Hz. The sampling rate for recording was 200 Hz. In January-February 2018, the 3-D seismic reflection experiment was conducted by the U.S. R/V Marcus G. Langseth along 58 seismic lines with an average spacing of 300 m in the same area of the OBS network (MGL1801 participants, 2018). For this survey, over 140,000 air gun signals were triggered every 25 m in flip-flop mode and covered a 60-km × 14-km area extending from the northern Tuaheni landslide complex and the easterly forearc slope to the subduction trench ( Figure 1c). The Langseth 36-element air gun array system was used as the seismic source. The air volume of each shot was 3,300 cubic inches and the towing depth was 7 m. Ninety-seven OBSs were successfully recovered by R/V Tangaroa in March-April 2018 (Kellett et al., 2018), all of which recorded the air gun signals generated from the Langseth.
After the recovery of the OBSs, the seafloor position of OBSs were determined using traveltimes of direct water waves. For locating OBSs, we used the method of Oshida et al. (2008) that comprises two steps of analysis. In the first step, a global search was performed at a grid interval of 0.0001°to find the node on the seafloor that best explains the observed arrival times of water waves. In the second step, the location determined by the first step was used as an initial point and a nonlinear inversion was carried out to improve the location estimation. The average 1σ location uncertainty calculated from traveltime misfits are~3 m. We then picked first arrivals manually. For picking, we used either hydrophone or vertical component data-whichever showed better data quality-and applied 3-to 12-Hz band-pass filter to improve the signal to noise ratio. The dense air gun shooting with 25-m interval designed for the 3-D reflection survey was not optimal for the OBS refraction study. Water waves from previous shots have high amplitude and overlap primary refraction phases from the next shot at horizontal distances between shot and receiver of~25 and~45 km, making it difficult to pick first arrivals at these offsets ( Figure 2). Otherwise, the OBS data are of good quality and P wave first arrivals can be traced up to 20-to 50-km distance for most OBSs. Since the air gun shooting was much denser than typical refraction studies and produced significantly redundancy in the data, we picked first arrivals from every fifth shot line for all stations. To improve the ray coverage at deeper parts, we additionally picked other shot lines for some stations that have a good signal to noise ratio and clearly recorded long-offset first arrivals (Figures 2 and S1 in the supporting information). We assigned picking uncertainty of 0.030 s to traces in the offset range of 0-15 km. This value of uncertainty was chosen by examining the wavelength of first arrivals observed in the offset range. Thereafter, we increased the uncertainty to 0.040 s between an offset of 15-20 km because many OBS records show a significant amplitude reduction in this offset range. The uncertainty was further increased to 0.050 s at 20-30 km offset because the apparent velocity of first arrivals increases to~6.0 km/s at this offset and the waveforms become more enriched in lower frequency compared to near-offset traces. Finally, the largest uncertainty of 0.060 s was assigned to traces >30 km offset where the first arrivals are possible to pick but contaminated by previous shot noise. For the tomography analysis, we used over 430,000 first arrivals.
In the seismic record sections, P waves with offset of <10 km arrive with apparent velocities less than 4.0 km/s, which indicates that the survey area is entirely covered by thick low-velocity sediments ( Figure 2). OBSs on the landward side of the trench show a significant reduction in amplitude of refraction phases at 10-to 20-km distances. This feature is commonly observed at the forearc region of subduction zones and may suggest that there are shadow zones due to velocity reversals and/or highly attenuative zones caused by complex forearc structure. At farther offsets, first arrivals are refractions with apparent velocities of~6.0 km/s for all the OBSs.

Tomography Analysis
We use 3-D azimuthal anisotropic inversion (Dunn et al., 2005) where the Vp model is parameterized with an isotropic component and two azimuthal anisotropy parameters. For each model grid, the anisotropic slowness, u ani , is defined as where u iso is the isotropic slowness and θ is the raypath azimuth. The terms A and B are anisotropy parameters. The magnitude of anisotropy, M, is expressed as or the percentage of anisotropy, α,

10.1029/2020JB020433
Journal of Geophysical Research: Solid Earth The fast direction, ψ, is The reduction velocity is 6.0 km/s. Although strong water arrivals from previous shots overlap the primary refraction phases, first arrivals are clearly observed in the horizontal distance of 0-20, 25-45, and greater than 50 km.
Theoretically, higher-order terms of even multiples of θ are also included in the formula (Backus, 1965). While previous studies show that 4θ variation is recognizable in traveltime data, it is difficult to constrain its spatial variation since it is roughly equivalent to or smaller than the traveltime pick uncertainty (Dunn, 2015). We also found a 4θ variation too small to model with our data, and we only examine the 2θ variation as described by the equation above.
We defined a model space of 70 km in the trench-normal direction (x axis), 16 km in the trench-parallel direction (y axis), and 18 km in the depth direction (z axis). The x axis in the trench-normal direction was oriented on seismic profile 05CM-04 . This model was discretized 0.2 km horizontally and 0.15 km vertically for forward calculation of traveltimes and raypaths. The upper surface of the model corresponds to the seafloor and was overlain by the water layer whose velocity was fixed at a constant value of 1.495 km/s.
Since the anisotropic velocity tomography is highly nonlinear and involves more parameters than isotropic analysis, the inverse calculation can be unstable even when a large quantity of data are available. To overcome this instability, we carefully constructed our starting model from existing, well-constrained velocity models and took multiple steps to sequentially develop toward anisotropic inversion. First, to obtain a stable starting model, a 2-D velocity model was created by taking the 1-D velocity profiles from two OBS refraction studies close to our study area, one from the forearc region and the other from the seaward side of the deformation front Kodaira et al., 2018) and linearly interpolating between them along the regional plate boundary model (Williams et al., 2013). The 2-D velocity model was then extended in the trench-parallel direction to create an initial 3-D model (Figure 3a).
No anisotropy was assumed in the initial model. Despite this careful preparation, the initial model has large data misfits (the root-mean-squared, RMS, residual of 558 ms and normalized chi-square misfit of 242.72), which may end up in a local minimum after inversions especially for the anisotropic parameters. To stabilize the calculation, we took the average of the resultant Vp model after three iterations and the Vp values of the starting model ("first" initial model) to make a smooth "second" initial isotropic model (Figure 3b). We repeated this process two more times (resulting in "third" and "fourth" initial isotropic models; Figures 3c and 3d) until the changes in the Vp values became insignificant. Although the inversion for the initial model preparation is anisotropic, the Vp results are similar to those from an isotropic inversion since we assigned weaker damping for Vp and stronger damping for anisotropy in these inversions.
At each round of inversion steps, we employed different node spacing for the inverse problem. For the first round, 2-3 times larger spacing were used than those used for the final round to first solve for larger-scale velocity variations. The node spacings used in the final round for the isotropic slowness was 0.5 km in the x axis and 0.6 km in the y axis. Vertically, the node spacing was 0.5 km at 0-to 8-km depths, 1.0 km at 8to 12-km depths, and 2.0 km from 12-km depth to the bottom of the model. The node spacing for anisotropic parameters were 1.5 km in the x and y axes and 0.5-2.0 km in the z axis. The inversion was regularized by smoothing and damping constraints that control the overall smoothness and stability of the final solution.
In this tomography code, the user needs to assign the strength of horizontal and vertical smoothness and the damping term on the magnitude of model changes for both isotropic slowness and anisotropy. We determined the values of all these inversion parameters by referring to our past experiences using the same code (Dunn et al., 2017) and by assessing the results for different inversion parameter tests (see supporting information).

Vp Structure
The final, preferred anisotropic model obtained after two iterations in the final round of inversion starting from the "fourth" initial model has the mean traveltime residual of 35 ms and the chi-square value of 0.86. Figures 4 and 5 show vertical cross sections and map view slices of the Vp values from this preferred model, respectively. One of the most important results is a low-velocity zone (Vp = 1.6-4.0 km/s) observed at the frontal part of the upper plate ( Figure 4). This zone forms an~30-km-wide and~5-km-thick accretionary wedge that is adjacent to a higher-velocity body (Vp = 2.0-4.5 km/s) on the landward side and corresponds to an area where multiple splay faults are developed ; Figure 1d). By contrast, no active faults are imaged in the higher-velocity region marking the seaward edge of the north Tuaheni basin. In map view, the high-velocity body is oblique to the trench, rapidly narrowing the width of the accretionary wedge from~30 km in the south to~20 km at the northern end of the 3-D box ( Figure 5).
Subducting along the Hikurangi Trough is the Hikurangi Plateau, a large igneous province of Cretaceous age with a crustal thickness over 10 km (Davy et al., 2008;Mochizuki et al., 2019). The structure of the incoming Hikurangi Plateau is different from typical oceanic crust and is characterized by thick low-velocity layers that shows a gradual increase in Vp with depth from~1.6 km/s at the seafloor to~5.0 km/s at 4-km depth below seafloor. This layer includes pelagic rocks of Paleogene to Miocene age and the Hikurangi basement mainly consisting of volcaniclastic rocks (Barnes et al., 2020). The lower volcaniclastic layer with Vp of 4.0-5.0 km/s extends in the downdip direction to beneath the accretionary wedge (Figure 4). Figure 6a shows plots of traveltime residuals calculated from the "first" initial model (Figure 3a), the final preferred model, and an isotropic model that was constructed in the same way as the preferred anisotropic model expect that only isotropic slowness parameters were inverted. This figure shows the best data fit for the preferred anisotropic model. Figure 6b indicates that many OBSs (those on the accretionary wedge and the north Tuaheni forearc basin in particular) exhibit clear variations in traveltime residuals that can be approximated by a sine function of 2θ in the isotropic model (green ticks in Figure 6b). These traveltime variations can be explained well by incorporating the anisotropic effect (red ticks in Figure 6b). The azimuthal variation in traveltime residuals becomes smaller on the seaward side of the trench (Figure 6b).    Figure 1c), on the accretionary wedge (indicated by white circles with a purple circumference in Figure 1c) and on the north Tuaheni forearc basin (indicated by white circles with a blue circumference in Figure 1c). The mean residuals and standard deviations were calculated for azimuth bins of every 20°and using the rays whose maximum depths are less than 6 km. The shot-receiver azimuth is measured from north. The X axis is oriented to the southeast-northwest direction (121°/−59°) as indicated by black arrows. While the isotropic model shows clear variations in traveltime residuals especially in the accretionary wedge and the forearc basin that can be expressed by a sine function of 2θ (green line), the preferred anisotropic model has better data fits independently of the azimuth direction.

Anisotropy Structure
Correspondingly, the 1-D depth profile of anisotropy shows that the accretionary wedge tends to show larger azimuthal anisotropy in the upper 2.0 km below the seafloor than the landward forearc basin and the seaward incoming plate (Figures 7 and S2). The map view slices reveal that areas with >5% anisotropy are found mostly in the accretionary wedge ( Figure 5). This 5% value is the average over the scale length our model can resolve (2-4 km); we cannot exclude the possibility that a much higher magnitude of anisotropy than 5% exists along the faults beyond our resolution. The cross section shows that the anisotropic features are not ubiquitous and homogeneous within the overriding plate and tend to become larger in the vicinity of the fault locations ( Figure 8).
The fast axes of Vp exhibit a clear spatial pattern: While they are mostly oriented in the trench-normal (i.e., plate convergence) direction in the accretionary wedge, they are rotated to the trench-parallel direction on the seaward side of the trench and in the landward north Tuaheni basin (Figures 5 and 9a). This regional variation is consistent with the results of shear wave splitting analysis (Zal et al., 2020; Figure 9b) and the regional stress field estimated by focal mechanism solutions (Townend et al., 2012). Occurrence of fast axes in the accretionary wedge agrees well with the directions of maximum horizontal stress inferred from borehole breakouts at Drilling Sites U1518 and U1519 (Figure 9a). These correlations demonstrate that our estimation of Vp anisotropy is consistent with other similar analyses and observations.

Model Resolution and Reliability
The reliability of obtained structural heterogeneities was evaluated by several different methods. Checkerboard tests using the same raypaths, matrix of partial derivatives, and smoothing/damping constraints as for the preferred model demonstrate that the anomaly pattern of 4 km horizontally and 3 km vertically is well recovered down to the depth of 6 km below seafloor ( Figure S3). Although slightly distorted and weaker at depths greater than 6 km, the checkerboard is recovered down to 9 km depth and we consider this to be the maximum depth of resolution in our velocity model. Smaller horizontal anomalies of 2 km are successfully reconstructed in the upper 4 km and partially recovered down to 6 km. These tests constrain the minimum size of velocity variations we can interpret. The checkerboard tests for anisotropy magnitude indicate that anisotropy is well resolved down to the depth of 4 km below seafloor and partially resolved down to 6 km beneath the accretionary wedge ( Figure S4). Because anisotropy is calculated from rays of different azimuths, our maximum resolution depth is limited by the maximum offsets recorded in the margin-parallel direction.
The robustness of our velocity model was also examined by checking the dependency on the choice of inversion parameters (see Figures S5-S10). The alternative smoothing parameters for isotropic slowness yield models that have only small velocity differences with respect to the preferred model (<0.1 km/s in most parts of the model) ( Figure S7). Different choices of smoothing parameters for anisotropy resulted in similar spatial pattern and magnitude in anisotropy (within 2% difference in most parts of the model), which strongly supports the validity of our anisotropy estimation ( Figure S8). We also performed similar tests for the damping strength of Vp and anisotropy and confirmed that the range of acceptable solutions is relatively small (Figures S9 and S10).
Finally, we compared our Vp model with the detailed Vp structure obtained by the FWI analysis of the 05CM-04 reflection profile  and logging data at two drill sites within the 3-D box   (Figure 10). Our velocity model has a spatial resolution of a few kilometers and thus is not sensitive to the abrupt velocity reductions in the upper~200 m of the model that are suggested by the FWI analysis and logging data, which have their vertical resolution of <~100 m. In addition, the

Journal of Geophysical Research: Solid Earth
shallowest part of our velocity model may be poorly constrained due to the lack of very near-offset first arrivals. Otherwise, our Vp model generally shows difference less than 0.2 km/s from the FWI and logging data at depths greater than 200 m, suggesting that the velocity variations in our Vp model are real within the range of overlap and likely good approximations deeper within the section.

Anisotropy Associated With Subduction Faults
Seismic anisotropy at shallow depths in subduction zones is an important indicator of the stress state and is commonly interpreted to be caused by preferentially oriented cracks in a compressional stress regime associated with plate subduction (Tonegawa et al., 2017;Zal et al., 2020). Our findings that (i) significant  (Figures S3 and S4). The coordinate system of the model (x and y axes) is shown in Figure 3.

10.1029/2020JB020433
Journal of Geophysical Research: Solid Earth anisotropy occurs associated with the active splay faults and the deformation front and (ii) the fast axes in the accretionary wedge lie parallel to the plate convergence direction suggest that the regional stress field and preferentially oriented cracks in the fault-bound blocks primarily control the seismic anisotropy. It is interesting that locations with large anisotropy tend to extend in the trench-parallel direction and some coincide with the steep slopes in seafloor bathymetry (Figure 9). Figure 7 compares the depth profile of anisotropy magnitude beneath the trench slope (that exhibits dip angles >20°) with that below the accretionary wedge with flat seafloor and shows that the former has slightly larger anisotropy in the upper 1 km. This difference may reflect a gravitational effect pointed out by Araragi et al. (2015) that a gravitational stress caused by topographic slopes influences the principal stress direction and thus the seismic anisotropy at shallow depths. Our results indicate that this effect may also play a role in trench slopes in subduction zones.
The maximum magnitude of anisotropy we have observed (>5%) is roughly equivalent to the values reported in oceanic spreading centers where the ridge-perpendicular tensile stress is dominant and thus likely to produce more cracks than at subduction zones (e.g., Dunn, 2015). The high-resolution velocity images produced by FWI revealed that some splay faults (but not all) were associated with low-velocity zones, suggesting that Figure 9. Map view slice of (a) anisotropy magnitude at 1-km depth below seafloor and comparison with (b) shear wave splitting anisotropy (Zal et al., 2020) and (c) seafloor depth. The fast axes of Vp agree well with the orientation of maximum horizontal stress (SH max ) inferred from borehole breakouts at Drilling Sites U1518 and U1519 green ticks) and shear wave anisotropy. At U1519, borehole breakouts with an approximate northeast-southwest (054°/234°) direction was observed between 600 and 650 m below seafloor . At U1518, borehole breakouts have an azimuth of 173°/353°from which a SH max orientation of east-is estimated. Blue ticks in (b) indicate the orientation of fast axes based on shear wave splitting (Zal et al., 2020). The seafloor bathymetry data in (c) was obtained during the 3-D seismic reflection cruise (Edwards et al., 2018). The study area includes the northern Tuaheni landslide complex (TLC north) and Tuaheni Ridge. A gray shade indicates the area of extension where normal faults across the upper slope of the north Tuaheni basin are distributed (Böttner et al., 2018). A black arrow in (c) shows the convergence direction at the Hikurangi trough between the Pacific Plate relative to the North Island forearc (Wallace et al., 2004). these faults may have high porosity and could act as fluid conduits . In the study area, it is expected from local seismicity that fluids released from the downgoing slab migrate into the overlying plate interface and upper plate (Shaddox & Schwartz, 2019;Todd et al., 2018;Warren-Smith et al., 2019). Such fluids contribute to increasing the crack density around the faults by keeping the cracks open at greater depths (e.g., Zatsepin & Crampin, 1997). However, not all splay faults are characterized by an across-fault velocity reduction and associated porosity increase. For example, while our Vp model indicates a high anisotropy area corresponding to the Pāpaku fault that is located in the trench slope ( Figure 11), a low-velocity volume is not associated with this particular fault in the FWI analysis of Gray et al. (2019). In addition, core samples from U1518 that crosscut the Pāpaku fault at~300-m depth suggest that there is not a significant increase in porosity at the fault zone . These considerations indicate that a factor other than preferentially oriented cracks may contribute to seismic anisotropy.
Although clay fabric is not observed in the core samples from the Pāpaku fault at U1518 , clay-rich sediments with evidence of ductile flow deformation are often observed at fault zones in subduction zones (Kinoshita et al., 2009;Tobin et al., 2001). In the northern Hikurangi margin, IODP drilling has confirmed that input sediments on the Hikurangi Plateau contain abundant clay minerals (Barnes et al., 2020). The sedimentary sequences are likely to stack and form coherent strata with a total thickness of up to several hundred meters in the vicinity of the fault zones during the process of accretion and/or subsequent deformation (Rowe et al., 2013). If the clay minerals in such thick sedimentary layers are preferentially aligned by the compressional stress by plate convergence, they may also result in significant upper plate anisotropy that can be detected by seismological methods.

Regional Stress Field Associated With Plate Subduction
Our velocity model revealed two other areas with different spatial orientation of anisotropy from the accretionary wedge which reflect variable stress field associated with plate subduction (Figure 5). First is the seaward side of the deformation front, where the fast axes of Vp are parallel to the strike of the trench at 2-to 3-km depths (Figures 5e and 5f). This feature is consistent with previous observation of anisotropy for the subducted Pacific Plate (Eberhart-Phillips & Reyners, 2009) and suggests that the maximum horizontal stress is oriented parallel to the trench. This stress regime is explained by the flexural bending of the incoming plate (e.g., Bodine et al., 1981). Although the overall magnitude is smaller (mostly less than 2%) compared to the magnitude of anisotropy in the accretionary wedge, this stress field may be active as implied by normal faults within the Hikurangi Plateau basement observed in seismic reflection data ( Figure 1d). On the other hand, the fast axes at 1-km depth is more variable. The drill core data at U1520 suggests that this depth corresponds to the uppermost part of the Cretaceous volcaniclastic layer consisting of the Hikurangi Plateau (Barnes et al., 2020) and seismic studies suggest that its top surface exhibits significant relief with horizontal scales of <10 km Barnes et al., 2020;Gray et al., 2019). We thus suggest that the scattered pattern of anisotropy in the shallowest part of the incoming plate may be a consequence of the stress field disturbed by structural heterogeneities at that depth.
The other area is the most landward part of the 3-D box where the fast axes of Vp rapidly change to the trench-parallel direction (Figures 5,9,and 12). This change may be related to the Tuaheni Ridge and the edge of the North Tuaheni basin, in an area where forearc extension has been previously mapped (Böttner et al., 2018;Figures 9c and 12a). The variation in seismic anisotropy is consistent with the results of shear wave splitting analysis (Zal et al., 2020; Figure 9b). While the depth variation in anisotropy was not well constrained by Zal et al. (2020), our results demonstrate that the upper few kilometers in the overriding plate is the region which determines this spatial pattern.
It is also intriguing that the spatial variation in anisotropy correlates well with the structural heterogeneities within the overriding plate. Figure 12 shows that in the north Tuaheni forearc basin significant anisotropy exists just above two features consisting of high-velocity material with 4.0-5.0 km/s. These bodies are likely to be the lowest part of the overriding plate and probably associated with the Mesozoic and Paleogene imbricated wedge backstop. Forearc normal faulting has been identified in the Nankai accretionary margins and suggested to result from decoupling of shallow strata from the underlying compressional wedge (Gulick et al., 2010;Moore et al., 2013;Sacks et al., 2013). Böttner et al. (2018) were unable to identify the cause of extension, but we find that spatial orientation in anisotropy correlates well with the structural heterogeneities within the upper plate and supports the suggestion of local uplift giving rise to the extensional strain.

. Implications on Seamount Subduction
On the basis of the seismic reflection and magnetic data, it has been suggested that an~20-km-wide seamount is subducting just beneath the central part of our 3-D survey area Bell et al., 2010; Figure 1). The location and characteristics of the seamount are important because many previous studies have discussed seismic activity in relation to the seamount (e.g., Todd et al., 2018). In addition, a subducting seamount may have a strong control on seismogenic behavior, hydrological processes and the structural evolution of the plate interface and the overriding plate (Sun et al., 2020). At the proposed location of the seamount, our velocity model imaged a low-velocity (Vp = 4.0-5.0 km/s) volume and did not detect any rapid velocity increase with depth ( Figure 11). This feature is consistent with the low seismic reflectivity at the top of the seamount while still accounting for the location of a high-relief magnetic source on the subducting plate .  (Böttner et al., 2018; a gray shade). We calculated the slope gradient from the bathymetric grid data using GMT (Wessel & Smith, 1998) Based on the currently available information, we propose two structural models focusing on the potential subducting seamount ( Figure 13). First, if the location of the seamount by Barker et al. (2018) is correct, our Vp model suggests that the internal structure of the seamount may be heterogeneous containing substantial regions with low seismic velocities (Figure 13a). Although this explanation seems apparently counterintuitive, because seamounts are usually topographic manifestations of igneous intrusions and thus exhibit higher velocities than the surrounding areas, subducting seamounts may be able to exhibit lower velocities in the upper few kilometers due to significant fracturing and/or comprising extrusive volcaniclastic lithologies (Hammer et al., 1994). Indeed, low-velocity seamounts with a similar size are reported in the Ryukyu subduction zone, southwest Japan, where the Amami Plateau subducts beneath the Eurasian plate  and along the Kermadec arc where the Louisville Ridge intersects the Tonga-Kermadec subduction zone (Bassett & Watts, 2015).
The alternative model is that the seamount is a higher-velocity body, located~10 km further downdip than previously thought (Figure 13b). This is the location where our velocity model shows high-velocity (Vp > 5.0 km/s) topography just below the plate interface. Correspondingly, the previous reflection studies revealed high-amplitude seismic reflectivity at the same place (Bell et al., 2010). The high reflectivity has been previously interpreted as low-velocity fluid-rich layers. Our velocity model may provide another interpretation for this reflectivity where top of the high-velocity seamount forms a thick reflective zone. This subduction model in Figure 13b suggests the seamount forms a low-velocity (Vp = 4.0-5.0 km/s) subducting channel in its wake that consists of sedimentary and volcaniclastic rocks in the process of subduction. However, this scenario would require a different interpretation of the magnetic data and the positive anomaly that is characteristic of nearby unsubducted seamounts . Future studies of subduction structure based on 3-D MCS reflection and OBS wide-angle reflection data, coupled with better constrained magnetic modeling, will provide more rigorous constraints of the subducting seamount.  Barker et al. (2018) is correct, our Vp model suggests that the seamount exhibits low velocities which would imply the seamount is fractured. (b) Alternatively, the seamount is a high-velocity body, located further downdip.

Conclusions
We quantitatively constrained the 3-D Vp values and azimuthal anisotropy in the northern Hikurangi margin using seismic refraction data obtained by the NZ3D OBS experiment. We found evidence for fault-bound anisotropy and discovered that large anisotropy is associated with active thrust faults along the plate interface and within the upper plate. The localized pattern and anisotropy of >5% may suggest that not only preferentially oriented cracks but also fault-bound clay-rich sedimentary layers contribute to producing the significant upper plate anisotropy. The fast direction of the Vp is consistent with the plate convergence direction in the accretionary wedge and rotates to a trench-parallel direction beneath the Tuaheni basin and seaward incoming plate. We also found velocity variations associated with the subduction of the Hikurangi Plateau. We suggest that the subducting seamount, proposed by previous studies, shows low velocities due to fracturing and/or extrusive volcaniclastic rocks. Alternatively, it is possible the seamount is located 10 km further downdip than previously thought.