Shallow Tectonic Stress Magnitudes at the Hikurangi Subduction Margin, New Zealand

Quantifying tectonic stress magnitudes is crucial in understanding crustal deformation processes, fault geomechanics, and variable plate interface slip behaviors in subduction zones. The Hikurangi Subduction Margin (HSM), New Zealand, is characterized by along‐strike variation in interface slip behavior, which may be linked to tectonic stress variations within the overriding plate. This study constrains in situ stress magnitudes of the shallow (<3 km) overriding plate of the HSM to better understand its tectonics and how they relate to larger scale subduction dynamics. Results reveal σ3: Sv ratios of 0.6–1 at depths above 650–700 m TVD and 0.92–1 below this depth interval along the entire HSM. Additionally, for depths below 650–700 m TVD, SHmax: Sv ratios of 0.95–1.81 in the central HSM and 0.95–2.3 in the southern HSM are estimated. These stress ratios suggest a prevalent thrust to strike‐slip (σ1 = SHmax) faulting regime across the central and southern HSM. In the central HSM, the presence of NNE‐NE striking reverse faults co‐existing with a modern σ1 (SHmax) aligned ENE‐WSW suggests that overtime the stress state here evolved from a contractional to a strike‐slip state, where the compressional direction changes from perpendicular (NW‐SE) to oblique (ENE‐WSW) to the Hikurangi margin. This temporal change in stress state may be explained by forearc rotation, likely combined with the development of upper plate overpressures. In the southern HSM, the modern WNW‐ESE/NW‐SE σ1 (SHmax) and pre‐existing NNE‐NE striking reverse faults indicate that stress state remains contractional and perpendicular (NW‐SE) to the Hikurangi margin overtime.


10.1029/2022GC010836
2 of 30 behavior of plate boundary faults, the origin and controls of diverse fault slip patterns, and to better assess seismic and tsunamigenic hazards along subduction zones (Huffman & Saffer, 2016;Riedel et al., 2016;Wu et al., 2019).
The characterization of stress magnitudes at the HSM is currently limited to relative stress magnitudes derived from earthquake focal mechanisms at seismogenic depths (Townend et al., 2012), and direct measurements of the minimum principal stress magnitudes (σ 3 ), vertical stress magnitudes (S v ), pore pressures (P p ) at shallow depths (<3 km) (Burgreen-Chan et al., 2016b;Darby & Ellis, 2001;Darby & Funnell, 2001), and stress regime in one borehole (Tuhara-1A) in the central HSM (HRT, 2000).Observations of relative stress magnitudes (≤60 km) by Townend et al. (2012) indicate a predominantly strike-slip and normal faulting regime along the HSM.P p measured from repeat formation tests (RFTs) and modular dynamic tests (MDTs), and inferred from drilling mud weights reveal shallow (<5 km) overpressures within the upper plate of the central HSM (Burgreen-Chan et al., 2016b;Darby & Funnell, 2001).High pore pressure in central and northern HSM are attributed to disequilibrium compaction of Miocene sediments and porosity reduction due to high horizontal compressive stresses associated with subduction of Hikurangi Plateau beneath the continental crust of North Island, New Zealand (Burgreen-Chan et al., 2016b;Darby & Funnell, 2001).σ 3 magnitudes determined from leak-off tests (LOTs) are less than or close to S v magnitudes (Burgreen-Chan et al., 2016b), suggesting variable normal, strike-slip, and reverse faulting regimes along the HSM.
In this study, we apply an indirect approach to constrain the three principal stress magnitudes along the shallow HSM crust using openly available borehole data.We discuss our findings in the context of understanding the upper plate tectonics within the HSM forearc.This study, in combination with stress orientation studies already completed for the HSM, provides a deeper insight into the variable tectonic behaviors associated with subduction margins and will serve as crucial information to assist in future hazard assessments of this region.

Geological Setting and Background
The HSM at the east coast of North Island, New Zealand, is a site of recent scientific investigation into the complexity of subduction dynamics.The western region of the HSM is composed of Triassic to Early Cretaceous basement rocks, referred to as the Torlesse Supergroup.These rocks formed as part of an ancient accretionary complex located on the southeastern margin of the Gondwana supercontinent (Bradshaw, 1989;Mortimer, 2004;Mortimer et al., 2014).Torlesse basement rocks consist mainly of indurated sandstone and mudstone (greywacke and argillite) with minor volcanics (Field et al., 1997).The basement is unconformably covered by Cretaceous deep-marine sedimentary rocks, represented by the now significantly deformed Matawai Group rocks, including the Koranga, Gentle Annie, Springhill, Karekare, and Glenburn Formations (Figure 2).
During the Late Cretaceous to Paleocene marine, fine-grained clastic sediments and limestones of the Tinui Group (including the Tangaruhe, Whangai, and Waipawa Formations) were deposited over the entire East Coast region (Figure 2) (Field et al., 1997).The Tinui group was formed during a time when the East Coast of New Zealand was a passive continental margin, following the cessation or slowing of Pacific-Gondwana subduction, the latter of which was likely due to the Hikurangi Plateau-Chatham Rise collision or the cessation of sea-floor spreading at the Osbourn Trough (Davy, 2014;Davy et al., 2008;Mortimer et al., 2019;Sutherland & Hollis, 2001).Above the Tinui Group, sediments are dominated by massive smectitic mudstones of the Late Paleocene to Eocene Wanstead Formation, and the Eocene to Oligocene Weber Formation marl (Figure 2).
The HSM developed in the early Miocene (∼25 Ma) by westward subduction of the oceanic crust of the Hikurangi Plateau beneath the continental crust of the North Island of New Zealand (Davy, 1992;Davy et al., 2008).During the Neogene the overriding Australian Plate underwent abrupt changes from fine-grained mudstones and limestones in the Paleogene to coarse grained facies such as thick conglomeratic sequences in the Miocene, accompanied by intense tectonic deformation of Late Cretaceous and Paleogene sediments (development of active thrust faults and folding), and the onset of andesitic volcanism in Northland (Ballance, 1976;Barnes et al., 2002;Field et al., 1997;Kamp & Furlong, 2006;Lamb, 2011;Mazengarb & Speden, 2000;Rait et al., 1991;Stern et al., 2006).
Three major tectonic phases have played a role in the development of the HSM from the Neogene to the present-day: 1.An early to mid-Miocene stage of contractional tectonics resulting in the development of extensive reverse faulting, folding, and tectonic uplift along the HSM (Bailleul et al., 2013;Barnes et al., 2002; Chanier  , 1999).The Miocene Tolaga Group rocks are deposited during this time and are dominated by the presence of bathyal mudstone and flysch, with units of paralic conglomerate and neritic sandstone and limestone also present (Figure 2) (Field et al., 1997).2. A mid-Miocene to early Pliocene stage of mixed contractional and extensional tectonics resulting in widespread contractional tectonic features across the HSM, and localized normal faulting and subsidence within the inner portion of the subduction wedge (Hawke's Bay region) (Barnes et al., 2002).The latter is attributed to tectonic erosion of the subduction margin (Chanier et al., 1999), or relief on the subducting Pacific plate (Barnes et al., 2002).In many areas, sedimentation appears to have continued uninterrupted across the Miocene-Pliocene boundary; however, early to mid-Pliocene erosion at some locations (e.g., Bland et al., 2008;Browne, 2004) has created unconformities that give the appearance of a sedimentation break across the Miocene-Pliocene boundary (Crutchley et al., 2016).3. A late-Neogene to present-day (Quaternary) stage (<10 Ma) of contractional tectonics associated with the structural inversion of listric extensional faults and folds, and rapid uplift of the Coastal Ranges, Axial Ranges, and Raukumara Peninsula (Bailleul et al., 2013;Barnes et al., 2002;Beanland & Haines, 1998;Nicol & Beavan, 2003;Nicol et al., 2002Nicol et al., , 2007)).The deformation in this stage is characterized by regional folding and mainly landward-dipping reverse faults (Bailleul et al., 2013).The Miocene to Pleistocene Mangaheia Group is characterized by rocks deposited in generally shallower-water environments than the underlying Tolaga Group.The Mangaheia Group rocks typically include shelfal mudstone, sandstone, and limestone facies, with some upper bathyal mudstone, and slope and basin floor fans (Figure 2) (Crutchley et al., 2016).
The oblique relative motion of the Australian-Pacific plate increases from 31 mm/year in the southern to 48 mm/year in the northern North Island (Figure 1a) (Wallace et al., 2004).The East Coast forearc has rotated at a rate of 3°-4°/Myr relative to the Australian plate, resulting in the TVZ back-arc rifting, strike-slip and/ or normal faulting in the onshore portion of the northern and central HSM, transpressional faulting in the southern HSM, and a large along-strike variation in convergence rate at the Hikurangi Trough (Figure 1a) (Fagereng & Ellis, 2009;Nicol et al., 2007;Wallace, Fagereng, & Ellis, 2012;Wallace et al., 2004).The oblique motion of the Australian-Pacific plate is partitioned into a margin-perpendicular component and a margin-parallel component.The margin-perpendicular component occurs along the Hikurangi subduction interface and provides NW-SE shortening mostly accommodated by slip on the subduction interface (>80%) and active frontal thrusts in the overriding plate (Nicol & Beavan, 2003).The margin-parallel component is largely accommodated by a combination of right-lateral strike-slip on the North Island Dextral Fault Belt (NIDFB) and clockwise rotation of the North Island forearc (Beanland & Haines, 1998;Nicol et al., 2007;Wallace et al., 2004).
Subduction of the buoyant Hikurangi Plateau (∼10-15 km thick) has been suggested as a potential cause of the uplift of a substantial part of the eastern North Island (the area between the volcanic arc and the trench) over the last few million years (Wood & Davy, 1994).This has resulted in the subduction plate boundary being located ∼40-120 km offshore and ∼12-15 km beneath the east coast of the North Island (Wallace, 2020;Williams et al., 2013).

Data Sources and Limits
Data used in this study are sourced from 44 boreholes along the HSM (Figure 3a), 41 of them are located within the onshore forearc and 3 are located offshore the east coast of NZ but west of the Hikurangi Trough.Data utilized includes wireline logging acquired over the period 1967-2013 from 0 to a maximum depth of 4,350 m below ground level.Wireline data include density logs from 27 boreholes (Table S1), sonic velocity logs from 25 boreholes (Table S1), and borehole image logs from 11 boreholes (Figure 3a).The positions of boreholes with borehole image log data relative to the Hikurangi subduction interface are shown in Figure 3b.Data presented here include the analysis of 18 LOTs and 39 formation integrity tests from 30 boreholes spanning a depth range of 71.5-3610.6 m (Table 1), mud weight logs from 44 boreholes, and RFT results from 2 boreholes spanning a depth range of 1,335-2,700 m.How each of these data are utilized in determining aspects of the in situ stress magnitudes across the HSM is detailed below.All depths in this study are referenced to ground level for onshore boreholes and sea level for offshore boreholes.

Vertical Stress Magnitude (S v )
Assuming that the vertical stress (S v ) is aligned to one of the principal stresses, the S v magnitude at any specific subsurface depth can be determined by the integration of rock densities from the surface to the depth of interest (Equation 1) (Bell, 1996;Engelder, 2014;King et al., 2010;Tingay et al., 2003): where ρ w is the average seawater density (1.03 g/cm 3 ),   is the gravitational acceleration constant (∼9.81 m/s 2 ), Z w is the depth of the water column (m), Z is the depth of interest (m), ρ(Z) (g/cm 3 ) is bulk density of the rock as a function of depth, and  ′  (g/cm 3 ) is the average density of the rock column above Z and below Z w .For onshore boreholes, Z w is equal to zero.
We utilize 26 density wireline logs to estimate S v profiles.At times wireline density logs are not acquired within the top depth intervals of drilled boreholes, the rock density is extrapolated from the top of a density log to the surface (seafloor for offshore boreholes) to more accurately determine a complete S v profile.This study uses several extrapolation methods: (a) using wireline sonic logs to convert compressional velocity to density values in boreholes where checkshot data or vertical seismic profile surveys are available (Kereru-1, Hawke Bay-1, Opoutama-1, Whakatu-1, Ngapaeruru-1, Tawatawa-1, and Titihaoa-1) (Bailey & Henson, 2019;King et al., 2010;Tingay et al., 2003), (b) using average densities from nearby boreholes with similar stratigraphy and tectonic setting (e.g., boreholes Kauhauroa-1, Kauhauroa-2, Kauhauroa-5, Makareao-1, and Tuhara-1A are all within <20 km of each other) (Traugott, 1997), or 3) using Gardner's relationship (Gardner et al., 1974) and/or regional Gardner's relationship (Table S2 in Supporting Information S1) to convert compressional velocity data from sonic wireline logs to density data logs (e.g., Hawke Bay-1, Rere-1) (Tingay et al., 2003).All density logs used in this study, supplied by the New Zealand Petroleum and Minerals group (NZPM), have undergone borehole environmental corrections.The density-logging tool needs contact with the borehole wall to yield accurate measurements (Asquith & Gibson, 1982).Therefore, rugosity of the borehole wall can cause poor tool contact with the formation, resulting in the density log measuring the density of the drilling mud and mud cake instead of the formation.Thus, it is necessary to correct density logs for spurious readings prior to calculating vertical stress (Bell, 1996).The density correction curve (DRHO), which reflects the difference between short-and long-spaced density measurements, is commonly used to assess the quality of bulk density data (Tingay et al., 2003).DRHO values greater than 0.2 g/cm 3 are typically considered as unreliable density data (Asquith & Gibson, 1982).Finally, the density log was manually de-spiked where sharp jumps occurred during cycle skipping (Tingay et al., 2003) and smoothed over 3 m interval to remove anomalous measurements.

Minimum Principal Stress Magnitude (σ 3 )
σ 3 can be measured directly from pressure-time plots produced during LOTs, extended LOTs (XLOTs), or minifrac tests (Addis et al., 1998;Bell, 2003;White et al., 2002;Zoback et al., 2003).In the HSM, LOTs are the most common tests available to calculate in situ  3 magnitudes.LOTs are pumping pressure tests conducted in a borehole a few meters below recently set casing shoes.During constant fluid volume pumping, the recorded fluid pressure increase stops behaving linearly with time as the injected fluid pressure surpasses the σ 3 confining stress around the borehole and fluid starts to penetrate into the formation around the borehole (Addis et al., 1998;Bell, 1996).The point when the fluid pressure-time curve becomes non-linear (leak-off pressure (LOP)) can be read as an approximation of σ 3 magnitude.If a LOT is stopped at any point before the LOP is reached, the test is called a FIT and fluid pressure has not exceeded σ 3 magnitude.In this case, the final fluid pressure value recorded during the FIT can be used as an estimate of the lower boundary of the σ 3 magnitude (e.g., Makareao-1, Zoback et al., 2003).
In the majority of boreholes studied here, the validity and accuracy of LOTs cannot be assessed as the pressuretime record data is not fully reported, with only the final LOP being provided in the text reports by drilling companies.Furthermore, pressure-time records are sometimes estimated by only a few distinct data points, obtained from pressure measurements on fluctuating gauges or flow rate estimations from counting pump strokes, making it impossible to determine the specific and accurate LOP values (Zoback, 2007).It is therefore possible for LOP to be reported slightly higher or extremely close to S v when the measurements are not carefully taken or reported.Further consideration for subduction margins is provided by Couzens-Schultz and Chan ( 2010), who demonstrate that in active compressional settings and seismically active regions, LOTs cause shear failure along pre-existing fractures rather than generating new tensile fractures, leading to an underestimation of the σ 3 magnitude.
We first calculate σ 3 :S v for all boreholes for which LOP measurements are available and then use the average of these data to extrapolate the σ 3 values beyond the depth of measurements.The σ 3 : S v = 1 is used to define the upper limit of the σ 3 profile in each borehole.

S Hmax Estimation From Borehole Failure Analysis
When a vertical borehole is drilled into a homogeneous, isotropic, and elastic medium parallel to one of the three principal stress orientations, the stress at the borehole wall is redistributed regarding to non-uniform, far-field principal stresses (Jaeger et al., 2009;Zoback, 2007).Assuming that the far field principal stresses are vertical and horizontal, the local principal effective stresses at a vertical borehole wall can be defined as follows (Moos & Zoback, 1990;Zoback, 2007): where σ θθ is the effective hoop stress (acting parallel to the borehole wall), σ ZZ is the effective vertical stress, σ rr is the effective radial stress (acting perpendicular to the borehole wall), S Hmax and S hmin are the maximum and minimum horizontal principal stress magnitudes, ϑ is Poisson's ratio, APRS is the annulus pressure at the time of borehole failure (or mud weight pressure), P p is pore pressure, and θ is the angle between the edge of borehole breakout and the S Hmax orientation (Figures 4a and 4b).
Where local effective stresses exceed the tensile or compressive rock strength of the formation around the borehole, borehole failures such as drilling induced tensile fractures (DITFs) and borehole breakouts (BOs) can form, respectively (Figure 4a).Measurements of the properties of these borehole failures, for example, the azimuth angle of BOs and/or DITFs and the angular width of BOs, can be used to determine in situ principal stress orientations and to calculate in situ stress magnitudes present at the time of drilling.
DITFs typically appear as narrow, conductive (on resistivity image logs) or low amplitude and slower travel time (on acoustic image logs) pairs, ∼180° from each other around the borehole wall circumference (Figure 4c).DITFs are generally parallel or slightly inclined to the borehole axis in vertical to semi-vertical boreholes (Brudy & Zoback, 1999;Zoback, 2007).Where DITFs are observed, the magnitude of the far-field S Hmax can be constrained using Equation 3 (Zoback, 2007): where S Hmax and S hmin are maximum and minimum horizontal principal stresses, respectively, T 0 is the formation tensile strength, P p is pore pressure, APRS is annulus pressure (or mud weight), and σ ∆T is thermal stress arising from the difference between the drilling mud temperature and formation temperature.σ ΔT is applied where there is a noticeable difference between mud and rock temperature, such as geothermal boreholes.The tensile rock strength in sedimentary rocks is often quite small (a few MPa) and can be assumed to be zero in the analysis of DITFs (Brudy & Zoback, 1999).In this study, σ ΔT is considered negligible.
BOs form as enlargements of the borehole diameter on opposite sides of the borehole wall where σ θθ is large enough to exceed the formations compressional strength (Figure 4a) (Bell & Gough, 1979;Zoback, 2007).The σ θθ magnitude reaches a maximum at θ = ±90° (Figure 4a), which occurs at a borehole azimuth oriented perpendicular to the S Hmax direction (Figure 4b).BOs typically appear as a pair of wide, out-of-focus, conductive (in water-based mud; Figure 4d) or resistive (in oil-based mud) zones on resistivity image logs, or as zones of low acoustic amplitude and slower travel time on acoustic image logs.BOs are located ∼180° from each other around the circumference of the borehole wall (Figures 4b and 4d).S Hmax magnitudes can be estimated by measuring BO widths (W bo ) from borehole image logs using Equation 4 (Barton et al., 1988;Vernik & Zoback, 1992): where W bo is the angular width of the BO, UCS is unconfined/uniaxial compressive strength of the formation, P p is pore pressure, APRS is annulus pressure or mud weight, S hmin is the minimum horizontal principal stress magnitude, and σ ΔT is the thermal stress effect resulting from the difference between the drilling mud temperature and formation temperature.In this study, σ ΔT is considered negligible.
UCS is a key parameter in estimating S Hmax magnitude (Equation 4) and can either be directly measured from laboratory strength tests on core samples or estimated using empirical relationships between UCS and other rock properties (Chang et al., 2006).Direct measurements of rock strength are rare for the HSM.Borehole Waingaromia-2 in the northern HSM is the only borehole where a laboratory strength test was conducted on calcareous claystone and mudstone core samples (acquired at 132 and 362 m measured depth, respectively), providing UCS values of 1.1-1.2MPa and friction angles of 20.5°-32.1°(friction coefficient 0.37-0.64)(Indo-Pacific Energy (NZ) Ltd., 2002).However, no relationship between P-wave slowness (Δt c ) and UCS was established because no geophysical logs were obtained and velocity measurements on core samples are unavailable.Therefore, in this study, UCS values are indirectly estimated using empirical relationships between rock strength and Δt c .Empirical equations have been developed for different rock types, relating various rock properties to UCS across the world.In this study we utilize a variety of empirical relationships between UCS and sonic velocity by matching appropriate equations to dominant lithologies encountered along each studied borehole in an effort to reduce uncertainty in UCS values and thus S Hmax magnitude values.Upper and lower bounds of the UCS are determined using various published empirical relationships (Chang et al., 2006;GMI, 2010;Horsrud, 2001;Lashkaripour & Dusseault, 1993;McNally, 1987) to provide a range of possible S Hmax magnitudes (Figure 5).Details on the equations used in individual boreholes to determine the lower and upper limits of UCS can be found in Tables S3,  S4, and S5 in Supporting Information S1.
A further important parameter required to calculate S Hmax magnitudes and effective stresses is P p .Direct P p measurement tests such as RFTs and MDTs are the most reliable measurements (Gunter & Moore, 1986;Zoback, 2007).However, these direct Pp measurements are difficult to acquire, particularly in low permeability formations, and are often only conducted at depths where possible overpressures may exist (Dutta et al., 2021;Lee et al., 2022;Y. Z. Ma & Holditch, 2015;Zoback, 2007).Drilling mud weight logs can provide indirect, continuous approximations of the Pp along a borehole, and can be used as a proxy of P p assuming the mud weights have been chosen to stabilize the borehole during drilling, and if no significant mud losses or kicks are reported (Van Ruth et al., 2002).Mud Losses of greater than 25 bbl/hr for water-based mud (Zhang & Yin, 2017) may indicate that annulus pressure exceeded Pp or/and σ 3 value, resulting in the loss of fluids into the formation.While kicks and high fluid influx indicate that Pp is greater than the annulus pressure.In both cases, the Pp derived from drilling mud weight logs should be corrected to generate a good estimation of Pp.In this study, we use mud weight logs from 41 boreholes to calculate Pp.Minor seepage (mud losses <22 bbl/hr) is reported for boreholes Kauhauroa-2/5, Makareao-1, Tuhara-1A, Ngapaeruru-1, Tawatawa-1, Titihaoa-1, and Ranui-2 in the intervals where BOs are observed, providing confidence in the use of mud weight logs in those intervals for Pp determination.A minor mud loss of 28 bbl/hr has been observed at severely fractured depth interval of 1,030-1,225 m TVD in borehole Ngapaeruru-1, which was easily treated by remedial techniques and procedures.Moreover, minor background gas and fluid influx are reported in boreholes Kahauuroa-5, Makareao-1, Tuhara-1A, Tawatawa-1, and Titihaoa-1, which were controlled by mud weight such that they never flowed.Since no significant mud losses or kicks are reported in the depth intervals where BOs are observed, we consider the annulus pressure records a good proxy of Pp in those depth intervals.The P p calculated from mud weight logs in Kauhauroa-5 and Titihaoa-1 boreholes are further calibrated using direct P p measurements obtained from RFTs.Formation tests conducted in 17 further HSM boreholes (Awatere-1, Hukarere-1, Kauhauroa-2/3/4B, Kiakia-1A, Makareao-1, Mangaone-1, Morere-1, Opoutama-1, Ruakituri-1, Takapau-1, Te Hoe-1, Tuhara-1 A/1B, Waitahora-1, and Waitaria-2) are not included in this study due to incomplete pressure build ups during testing in low-permeability formations, test seal failures, or tests conducted in formation intervals supercharged to hydrostatic pressure.

S Hmax Magnitude Estimation From Frictional Limit Theory
To constrain S Hmax magnitudes that result in the observed BO and DITF occurrences, the stress state is assumed to be limited by Coulomb frictional sliding on an optimally oriented and pre-existing fault plane (Zoback, 2007).This means that the maximum effective principal stress cannot exceed the stress value required to cause slip, defined by the friction coefficient (μ) of adjacent faults, on a critically oriented fault plane (Jaeger et al., 2009;Sibson, 1974): where σ 1 is the maximum principal stress, σ 3 minimum principal stress, P p is the pore pressure, and μ is the coefficient of friction on an optimally oriented, cohesionless, pre-existing fault.
Recent stress magnitude studies recommend using modified Wiebols-Cook, modified Lade criteria, or modified Drucker-Prager failure criteria for sand and shale layers, a brittle/frictional failure criterion that incorporates the effects of the intermediate principal stress on rock strength.However, these methods require polyaxial rock strength measurements and intermediate principal stress measurements (Chang et al., 2010;Chang & Song, 2016;Huffman & Saffer, 2016;W Lin, 2014), the former of which are not available for the region explored in this study and the latter of which are not well constrained.
This constraint is typically displayed as a stress polygon, which shows the permissible values of horizontal principal stress magnitudes for a specific depth, S v , μ, and P p for normal, strike-slip, and thrust faulting tectonics (Zoback, 2007).Although this method only provides the upper and lower limits for the S Hmax magnitude, it can yield more accurate ranges of permissible S Hmax magnitudes when combined with S Hmax magnitude estimates from borehole failure analysis (Chang et al., 2010;Huffman & Saffer, 2016).

Tectonic Stress Regime Index (A 𝜙 )
In order to characterize a stress regime or faulting style with stress magnitude data, we use the stress regime index (A  , Equations 6a and 6b) described by Simpson (Delvaux et al., 1997): where n is the number of principal stress components greater than the principal stress whose axis is closest to the vertical, R is the stress ratio, and σ 1 , σ 2 , σ 3 are the maximum, medium, and minimum principal stress magnitudes, respectively.
A  values range from 0 to 1 in normal faulting regimes, 1 to 2 in strike-slip regimes, and 2 to 3 in thrust faulting regimes.

Vertical Stress Magnitudes
S v magnitudes determined from 24 onshore boreholes provide overburden stress gradients ranging from 20.92 to 26.97 MPa/km, with a mean value of 22.58 ± 1.23 MPa/km (Figure 6a; Table S2 in Supporting Information S1). S v magnitudes measured within the 3 offshore boreholes range from 20.9 to 21.7 MPa/km (Figure 6b) with a mean value of 21.26 ± 0.4 MPa/km.

Minimum Principal Stress Magnitudes
σ 3 magnitudes calculated from 21 LOT data range from 1.9 MPa at 102.3 m TVD (borehole Waitaria-2) to 77.9 MPa at 3610.6 m TVD (borehole Rere-1) (Table 1).For all examined LOT data (with the exception of one test in borehole Titihoa-1 where σ 3 derived from LOP is greater than S v ), the normalized σ 3 ratio ranges from 0.6 to 1 (Table 1).It is not surprising to find that LOP is greater than S v , as it is commonly suggested by studies that LOP generally tends to slightly overestimate σ 3 (Addis et al., 1998;Bell, 1996;Breckels & van Eekelen, 1982;White et al., 2002).Therefore, in this study, we correct LOP > S v to LOP = S v .
LOP values measured in boreholes Tuhara-1/1A and Ranui-1 (11 m west of Ranui-2) are remarkably low, such that LOP values are less than the σ 3 values estimated by normal faulting failure with friction coefficients less than 0.6.In borehole Titihoa-1, the σ 3 value (43.4 MPa) derived from an LOT performed at 1979.8 m TVD is greater than the S V for this depth (41.3 MPa).In this case, σ 3 is considered to be vertical, indicating a thrust/reverse faulting regime (Zoback, 2007).In borehole Tawatawa-1, two LOTs were performed at 722.5 m TVD.The initial test yielded a LOP of 13.23 MPa, while the second LOT yielded a LOP of 13.36 MPa (Tap Oil Limited, 2004).Our reassessment of pressure-time curve of the second LOT (which had more data defining the time-pressure plot) reveals that the formation breakdown pressure was reported rather than LOP, resulting in an overestimation of σ 3 making it appear greater than the S v for this depth.We determine the LOP of the second test by intersecting the straight line of the linear section with the tangent line of the ascending section on the pressure-volume curve (Figure 7) and report a σ 3 magnitude of 12.7 MPa, almost equal to S v (13 MPa).
Analysis of 39 FIT measurements show that FIT data are typically below S v in most boreholes in this study.However, in some boreholes (Table 1), FIT results are approximately equal to or greater than S v , suggesting σ 3 = S v .
We calculated σ 3 : S v (after correcting LOP > S v to LOP = S v ) for all 21 LOT measurements in the study area and found that σ 3 : S v varies significantly above and below 650-700 m TVD (Figure 8).The σ 3 : S v ranges from 0.6 to 1 at depths above 650-700 m TVD, while it ranges from 0.92 to 1 below this depth interval (Figure 8).Above and below 650-700 m TVD, the average σ 3 : Sv values derived from LOT measurements are 0.79 and 0.95, respectively.
In this study, to create the σ 3 profile in boreholes where LOT measurements are not available at the depth of interest, we only consider LOT and FIT = S V measurements recorded at depths below 650-700 m TVD for two main reasons: (a) to exclude the influence of topographic effects and shallow processes such as gravitational collapse, erosion, and subsidence in the σ 3 calculation, and (b) because our borehole breakouts data to estimate S Hmax magnitudes in 7 out of 8 boreholes are located below 700 m TVD.

Pore Pressure Magnitudes
In this study, P p is determined indirectly using drilling mud weights from 41 boreholes (all of which are drilled <5 km depth) and from direct measurements such as RFTs and MDTs obtained from two boreholes (Kauhauroa-5 and Titihaoa-1).The overpressure ratio (λ = Pp:S v ) ranges between 0.45 and 1 in the central HSM and between  0.45 and 0.85 in the southern HSM (Figure 9).The λ = 0.45 indicates a pore pressure equal to hydrostatic pressure (P hyd ), and the λ = 1 indicates a pore pressure equal to S v .Overpressure zones are present along the entire HSM and vary with respect to the depths they occur at (in part related to stratigraphy) and their magnitudes (Figure 9) (Burgreen-Chan et al., 2016a;Darby & Ellis, 2001).HSM overpressure zones are typically found in Cretaceous-Paleogene sequences overlain by the smectite-rich seal rocks of the Eocene Wanstead Formation, are spatially variable in Neogene units, and are rarely found in Quaternary sequences (Burgreen-Chan et al., 2016a;Darby & Ellis, 2001).Available data suggest that the northern and central HSM show larger overpressures at shallower depths compared to the southern HSM.Extremely high overpressure zones (λ ≥ 0.9) are observed in boreholes Rotokautuku-1, Kauhauroa-3/4 in the northern and central HSM at depths as shallow as 330 m below the ground level (Figure 9).

Stress Magnitude From Borehole Data
We analyze borehole image logs acquired from 11 boreholes; however, BOs and DITFs are only observed in 10 of these image logs (Figure 1b).Among the remaining 10 boreholes, two boreholes (boreholes 8 and 10; Figure 1b) are only drilled to depths <200 m from the surface; therefore, we exclude them from this study as stress features forming in these boreholes may be related to more complex forces beyond that of the horizontal principle stresses.A range of potential S Hmax magnitudes are estimated along 8 boreholes at depth intervals where BO widths and DITFs are measured from their image logs.Lower and upper limits of S Hmax magnitudes are constrained by theoretical limits provided by slip on pre-existing faults with a friction coefficient of 0.6 (described in Section 3.4.2).

Makareao-1 Borehole
The FIT: S v ratios of 1.03 and 1 are determined at depths 306.2 and 484.8 m TVD, respectively (Table 1).The FIT: S v ≥ 1 indicates that σ 3 = S v in this borehole.The S Hmax :S v ratio of 1.03-1.31 is determined for borehole Makareao-1 using the S Hmax values calculated from the lower and upper values of UCS.The S Hmax :S v ratio of 1.03-1.31and a σ 3 : S v = 1 along this borehole indicates a stress regime such that σ 3 = S v < S Hmax (Figure 10a).A ϕ = 2 is determined from calculated stress magnitude data in this borehole.

Kauhauroa-2 Borehole
A σ 3 : S v ratio of 0.81 is determined from σ 3 value calculated using LOT data at 463.8 m TVD (Table 1).The σ 3 values in the deeper part of the borehole are calculated from the average σ 3 : S v ratio of 0.95 along the HSM and are further constrained by the lower limit of σ 3 value determined from an FIT = 30.23MPa at 1707.3 m TVD (Table 1).The S Hmax : S v ratio of 0.95-1.71 is determined for borehole Kauhauroa-2 using the S Hmax values calculated from the lower and upper value of UCS (Figure 10b).The S Hmax : S v ratio of 0.95-1.71and the σ 3 : S v ratio of 0.95 indicate a dominant stress regime such that S hmin ≤ S v ≤ S Hmax (Figure 10b).A 0 ≤ A ϕ ≤ 1.94 is determined from calculated stress magnitude data in this borehole.S v , S hmin , and the lower limit of S Hmax are nearly equal below 1980 m TVD such that S hmin ≈ S Hmax ≈ S v .

Kauhauroa-5 Borehole
A σ 3 :S v ratio of 0.73 is determined from σ 3 value calculated using LOT data at 459.2 m TVD (Table 1).The σ 3 values in the deeper part of the borehole are calculated from the average σ 3 :S v ratio of 0.95 along the HSM and are further constrained by the lower limit of σ 3 value calculated from an FIT value of 27.13 MPa at 1276.1 m TVD.The S Hmax :S v ratios of 0.95-1.13 in 1280-1350 m TVD and 0.97-1.54 in 1390-1750 m TVD are determined in this borehole using the S Hmax values calculated from the lower and upper values of UCS.
The analysis of S Hmax magnitudes and the σ 3 :S v ratio of 0.95 indicate a dominant S hmin ≈ S v ≈ S Hmax (S Hmax-S hmin < 5 MPa) and 0 ≤ A ϕ ≤ 1.13 in the depth interval of 1280-1350 m TVD (Figure 10c).Moving along the depth to 1390-1750 m TVD, S hmin ≤ S v ≤ S Hmax and 0.44 ≤ A ϕ ≤ 1.92 are observed.
Further constraints on stress magnitudes are made in this borehole using the presence of DITFs between 1741 and 1745 m on FMI borehole image logs (Figure 10c).The presence of DITFs at 1742.3 m suggests that the S Hmax should be above the DITF line (Figure 11), where the local hoop stress can be tensile (Equation 3), but also inside the stress polygon with μ = 0.6.The possible range of S Hmax constrained using this information lies inside the blue shaded area (Figure 11), suggesting a S Hmax = 49.2-61MPa and a stress state such that S hmin ≤ S v < S Hmax .

Tuhara-1A Borehole
A σ 3 : S v ratio of 0.63 is determined from σ 3 value calculated using LOT data at 590.8 m TVD (Table 1).The σ 3 values in the deeper part of the borehole are calculated from the average σ 3 :S v ratio of 0.95 along the HSM, and are further constrained by the lower limit of σ 3 value determined from the FIT value of 40.6 MPa at 2149.5 m TVD.The S Hmax :S v ratio of 0.95-1.81 is determined for borehole Tuhara-1A using the S Hmax values calculated from the lower and upper values of UCS.
The S Hmax :S v ratio of 0.95-1.81and the σ 3 :S v ratio of 0.95 indicate a dominant stress regime such that S hmin ≤ S v ≤ S Hmax (Figure 5).A 0 ≤ A ϕ ≤ 1.95 is determined from calculated stress magnitude data in this borehole.

Ngapaeruru-1 Borehole
The σ 3 values in this borehole are calculated from the average σ 3 :S v ratio of 0.95 along the HSM and are further constrained by the lower limit of σ 3 value determined from FIT values of 8.35 and 16.86 MPa at 501.9 and 962.7 m TVD, respectively.The S Hmax :S v ratio of 0.95-1.75 is determined for borehole Ngapaeruru-1 using the S Hmax values calculated from the lower and upper values of UCS.
The S Hmax :S v ratio of 0.95-1.75 and the σ 3 :S v ratio of 0.95 indicate a dominant stress regime such that S hmin ≤ S v ≤ S Hmax (Figure 12a).A 0 ≤ A ϕ ≤ 1.94 is determined from calculated stress magnitude data in this borehole (Figure 12a).The upper limit of S Hmax magnitudes from the upper values of UCS are constrained by the limits provided by slip on pre-existing faults with μ = 0.6.

Tawatawa-1 Borehole
A σ 3 :S v ratio of 1 is determined from σ 3 value calculated using LOT data at 722.5 m TVD (Table 1).The S Hmax :S v ratio of 1-1.82 is determined for borehole Tawatawa-1 using the S Hmax values calculated from the lower and upper values of UCS.The S Hmax :S v ratio of 1-1.82 and the σ 3 :S v = 1 indicate a dominant stress regime such that σ 3 = S v ≤ S Hmax (Figure 12d).The upper limit of S Hmax magnitudes from the upper values of UCS are constrained by the limits provided by slip on pre-existing faults with μ = 0.6.A ϕ = 2 is determined from calculated stress magnitude data in this borehole.

Titihaoa-1 Borehole
The σ 3 :S v ratios of 0.86, 0.94, and 1.05 are determined from σ 3 values calculated using LOT data at 614, 1585.7, and 1979.8 m TVD in this borehole (Table 1).The σ 3 :S v ratio of 1.05 at 1979.8 m TVD indicates that σ 3 = S v at depth intervals of 2,200-2,700 m TVD.The S Hmax :S v ratio of 1.02-1.41 is determined for borehole Titihaoa-1 using the S Hmax values calculated from the lower and upper values of UCS.
The analysis of S Hmax magnitudes and σ 3 :S v = 1 at depth intervals of 2200-2700 m TVD indicate a stress regime such that σ 3 = S v ≤ S Hmax (Figure 12c).The upper limit of S Hmax magnitudes from the upper values of UCS are constrained by the limits provided by slip on pre-existing faults with μ = 0.6.A ϕ = 2 is determined from calculated stress magnitude data in this borehole.

Ranui-2 Borehole
The σ 3 profile in this borehole is calculated from the average HSM σ 3 :S v ratio of 0.95 are further constrained by the lower limit of σ 3 value determined from FIT value of 6.35 MPa at 395 m TVD.The S Hmax :S v ratio of 0.95-2.3 is determined for borehole Ranui-2 using the S Hmax values calculated from the lower and upper values of UCS.
The S Hmax :S v ratio of 0.95-2.3and the σ 3 :S v ratio of 0.95 indicate a dominant stress regime such that S hmin ≤ S v ≤ S Hmax (Figure 12b).A 0 ≤ A ϕ ≤ 1.96 is determined from calculated stress magnitude data in this borehole.

Shallow HSM Tectonics
Stress magnitudes calculated from borehole data indicate that the S Hmax :S v ratios ranging from 0.95 to 1.81 in the central HSM and 0.95-2.3 in the the southern HSM.Additionally, σ 3 :S v ratios of 0.6-1 are estimated at depths above 650-700 m TVD, while 0.92-1 are estimated below this depth interval along the HSM.These stress magnitude results reveal that across the central and southern HSM, S Hmax is dominantly σ 1 , indicating a thrust to strike-slip faulting regime.The observed dominant thrust to strike-slip faulting regime is consistent with observed contractional tectonics in the HSM developed by the subduction of the Hikurangi Plateau beneath the North Island (Barnes et al., 1998;Nicol & Beavan, 2003), and the strike-slip faulting generated by forearc rotation of the East Coast (Beanland & Haines, 1998;Litchfield et al., 2014;Nicol et al., 2007;Wallace et al., 2004).Behboudi et al. (2022) report a dominant ENE-WSW shallow crust S Hmax orientation within the central HSM, and WNW-ESE to NW-SE S Hmax orientations for the southern HSM (Figure 1b).Considering σ 1 = S Hmax along the HSM, observed S Hmax orientations suggest that the contemporary maximum compressional stress switches from oblique (ENE-WSW) to the Hikurangi margin in the north and central HSM to roughly perpendicular (WNW-ESE to NW-SE) to the Hikurangi margin in the southern HSM.Based on our confirmation here that σ 1 = S Hmax along the HSM, it is likely that contemporary tectonics in the central HSM are dominantly strike-slip, while in the southern HSM, more contractional tectonics may be expected.
The NNE/NE striking faults in the central HSM, while currently inactive, express reverse dip-slip components to them based on seismic survey data (Western Energy New Zealand, 2001).These contractional faults developed in the early Miocene (∼25 Ma) as a result of the westward subduction of the Hikurangi Plateau beneath the continental crust of the North Island of New Zealand (Davy, 1992;Davy et al., 2008).This tectonic slip is at odds with the contemporary fault strike-parallel σ 1 (S Hmax ).These central HSM faults formed in an initially contractional stress state such that σ 3 = S v , σ 1 = S Hmax oriented NW-SE which would have been consistent with the NW-SE component of Pacific-Australian plate motion.We suggest that over geological time, this stress state changed from this contractional state to the modern strike-slip/contractional-oblique stress state (σ 3 : S v = 0.92-1, σ 1 = S Hmax oriented ENE-WSW).This switch in σ 1 orientation overtime and along HSM strike may be explained by (a) long-term clockwise rotation of the Hikurangi forearc (b) clockwise rotation of the Hikurangi forearc in conjunction with high shallow crust overpressures and/or mechanical property variations, and/or (c) along-strike variation in slip behavior in the HSM.
Clockwise rotation of the forearc, which accommodates the margin-parallel component of oblique Pacific-Australian plate motion, drives strike-slip and/or normal faulting within the onshore portion of the northern and central HSM, and transpressional faulting in the southern HSM (Figure 5, Fagereng & Ellis, 2009;Nicol et al., 2007;Wallace et al., 2004;Wallace, Fagereng, & Ellis, 2012).Behboudi et al. (2022) suggest that this forearc rotation is likely responsible for generating strike-slip stress state with ENE-WSW S Hmax = σ 1 in the central HSM, and contemporary contractional stress state with WNW-ESE/NW-SE S Hmax = σ 1 in the southern HSM.However, our stress magnitude results of σ 3 :S v = 0.92-1 and σ 1 = S Hmax leave a possibility for both strikeslip and contractional stress states to occur across both the central and southern HSM due to poorly constrained UCS values used in this study, a limitation of the study that could be restricted by laboratory rock strength testing of both onshore and offshore HSM lithologies.
The northern and central HSM have high Pp based on borehole data (Burgreen-Chan et al., 2016b;Darby & Funnell, 2001), magnetotellurics (Heise et al., 2019), and seismic tomography (Bassett et al., 2014;Eberhart-Phillips et al., 2017).Overpressure reduces the effective normal stress on fault planes, meaning that the existing NNE/NE striking faults in this region will be able to slip at lower shear stresses.Therefore, as a result of this overpressure, these faults could be less stable, allowing the hangingwall of upper plate faults to move more easily in response to NE-SW forces raised from forearc rotation.In this scenario, forces raised from forearc rotation were able to alter stress state overtime from σ 3 : S v = 1 and σ 1 = S Hmax with NW-SE S Hmax orientation, compatible with NW-SE component of Pacific-Australian plate motion and old geological structures, to σ 3 : S v = 0.92-1 and σ 1 = S Hmax with ENE-WSW S Hmax orientation.Similar shallow high overpressures are not observed in the hanging walls of upper plate faults in the onshore of the southern HSM.Therefore, it is possible that the NE-SW forces resulted from forearc rotation alone are insufficient to exceed the fault shear resistance and change the orientation of σ 1 away from the NW-SE component of Pacific-Australian plate motion, however they may have been high enough to play a role in reducing σ 3 magnitudes to the point that they become S v , resulting in a more transtensional tectonic regime overtime.
The mechanical properties of fault gauges and formations hosting faults (friction coefficient and rock strength) can play a role in controlling upper plate tectonic stresses (Mantovani et al., 2000;Marotta et al., 2002).Reiter (2021) investigated the impact of physical and elastic parameter contrasts on S Hmax orientation and proposed that contrasts in Young's modulus can introduce S Hmax rotations up to 78°, with larger stress rotations occurring within the softer lithologies.Behboudi et al. (2022) proposed that basement uplift in the southern HSM may introduce lateral geomechanical heterogeneities and variations in rock and sediment physical properties along the HSM, which may influence S Hmax orientations.Such that clay and sand-siltstone sediments (Miocene to present), where our stress data are calculated, in the upper plate of the central HSM may geomechanically differ from clay and sand-siltstone sediments (Miocene to Cretaceous) in the onshore of the southern HSM.Therefore, we analyzed physical properties of aforementioned sediments and discovered that P-wave velocity (V p ), UCS ranges, and dynamic Young's modulus (calculated from V s (S-wave velocity), V p , and density) in the central HSM are lower than the onshore of southern HSM (Figure 13).In this scenario, S Hmax orientations in the sediments of the central HSM, which have lower UCS and dynamic Young's modulus, could be easily reoriented in response to longterm forces such as forearc rotation compared to southern HSM.This theory, however, does not explain why the offshore boreholes in the southern HSM have not reoriented in response to long-term forearc rotation forces, while having comparable V p , UCS range, and dynamic Young's Modulus to the boreholes in the central HSM.
This along-strike variation in contemporary stress state is spatially consistent with north to south variation in slip behavior along the Hikurangi subduction interface (Figure 1a).In the northern and central HSM, the subduction interface is largely creeping and experiences shallow (<15 km), episodic SSEs that extend offshore and possibly to the trench.At the southern HSM, the plate interface is strongly interseismically locked to ∼30 km depth, and is currently accumulating elastic strain in the surrounding crust (Wallace, 2020).Some studies suggest that SSEs can release the amount of energy equivalent to a M w 6.5-8 earthquake (Dixon et al., 2014;Wallace, Beavan, et al., 2012).In the central HSM, the recurring SSEs and frequent earthquakes may release energy overtime such that the normal to shear stress ratio on pre-existing faults has changed in a way that makes it easier to slip in response to forces deriving from long-term forearc rotation.While stress accumulation due to locked nature of the southern HSM does not allow the normal to shear stress ratio to change considerably on the existing NNE/ NE striking compressional faults and makes it difficult for the hanging wall of these faults to slip in response to forearc rotation forces; therefore, stress state has not changed overtime in the southern HSM.However, the static stress drop of SSEs is estimated to range 0.01-1.0MPa (Gao et al., 2012).Given that the contemporary σ 3 : S v ≈ 0.95 and S v −σ 3 ranges between 0 and 3 MPa (for depths less than 3 km), these SSEs should have existed in the central HSM for more than 20 years such that they were able to release energy in order of 3 MPa (for depths less than 3 km) to change the initial σ 3 : S v = 1 to the contemporary σ 3 : S v = 0.92-1 and reorient the S Hmax orientation from NW-SE to ENE-WSW in this region.Further research and modeling are required to determine and quantify the initial stress state and whether the amount of stress released during SSEs in the central HSM was sufficient to support such a theory.

Extensional Tectonics Within the HSM Forearc
There are locales in the central and southern HSM where stress magnitude determination suggests a normal faulting regime (σ 3 :S v < 1 and 0 ≤ A ϕ ≤ 1).Also σ 3 :S v < 1 where σ 3 calculated from LOT data is observed for 13 tests conducted at depth intervals anywhere from ≈102 to 3611 m TVD in northern and central HSM boreholes (Table 1), and from 3 tests at depth intervals of ≈357-1,586 m TVD in southern HSM boreholes (Table 1).Several factors can result in localized normal faulting regime at subduction margins, including (a) uncertainties in calculated UCS values and/or σ 3 magnitudes used to determine stress states in this study, (b) the presence of local, active normal faults, and (c) fluctuations in stress magnitudes modulated by seismic cycles.

Uncertainties in Calculated Parameters for Stress Quantification
Estimations of S Hmax magnitudes are highly sensitive to the UCS values used, particularly when UCS is determined from empirical relationships not constrained by laboratory testing (Zoback, 2007).Due to the lack of direct UCS data in this region, and a lack of empirical relationships for the formations of this region to determine UCS from other rock properties, this study relied on the use or a range of empirical relationships developed elsewhere to generate a low and high limit for UCS at the HSM.These UCS ranges were then used to generate the lowest and highest limits of S Hmax magnitude.When the lowest limit of UCS is used, it can result in a potentially extensional stress state such that S hmin ≈ S Hmax ≤ S v .As such, the uncertainty in calculated UCS values, and the resulting potential errors it can introduce into a stress model for the HSM, highlight the importance of dedicated laboratory tests for developing robust empirical relationships for UCS in the HSM region, and subduction regions like this, where stress is a critical geological consideration for hazard and resource management.
Inaccuracy involved in LOT measurements (Section 3.3) along with the lack of detail reported on LOT results introduces an unknown level of uncertainty on estimated σ 3 magnitudes, and hence on estimated σ 3 :S v ratios using this data.Additionally, the lack of LOT data along each borehole necessitates the estimation of σ 3 profiles from the average σ 3 :S v = 0.95, which also carries uncertainty.As a result, we recognize the potential impact this has on calculations of S Hmax magnitudes here, as well as on any interpretations of regional stress state and tectonics.To investigate the potential effect of σ 3 uncertainties on S Hmax calculations, we use both the lower and upper limits of σ 3 values calculated from σ 3 :S v = 0.92-1, BO widths, and the lower and upper boundary of UCS values.This analysis reveals that the σ 3 magnitude uncertainties at the scale explored here have little influence on S Hmax magnitude calculations (±3.5 MPa) and hence do not change our findings about the stress regime and tectonics within the HSM (blue areas in Figure 5; Figures 10b and 10c; Figures 12a and 12b).
The σ 3 magnitude of 8.4 MPa measured from LOP in borehole Tuhara-1/1A (590.8 m TVD; Table 1) is lower than σ 3 values of 8.95 MPa estimated from normal faulting failure with a friction coefficient of 0.6 (Equation 5).This lower σ 3 magnitude may indicate that there are active normal faults at this depth along this borehole.In addition, borehole Tuhara-1A is located within the Tuhara anticline structure, formed by contractional stresses resulting from two blind thrust faults beneath the structure (Western Energy New Zealand, 1999).Our stress magnitudes and HRT's (2000) analysis suggest that the Tuhara structure currently experiences a dominant strike-slip faulting regime (S hmin ≤ S v ≤ S Hmax ; 1 ≤ A ϕ ≤ 2) along the majority of the borehole, interspersed with intervals of normal faulting regime (S hmin ≤ S Hmax ≤ S v ; 0 ≤ A ϕ < 1) mainly within the 1700-1820 m and 2100-2145 m TVD depth interval (Figure 5).A prominent feature of the Tuhara structure, as indicated by seismic reflection profiles, is the observation of relatively short steep east-and west-dipping normal faults throughout Pliocene and Miocene successions (Barnes et al., 2002;Western Energy New Zealand, 1999).Accordingly, we relate the appearance of normal stress states in our data to the normal structures that develop as part of the larger compressional structural architecture of this borehole site, and not due to the previously discussed uncertainties in the calculated UCS and/ or σ 3 magnitude values.This could particularly be the case where both the calculated lower and upper limits of S Hmax magnitudes are less than S V (e.g., at 1700-1820 mTVD in Tuhara-1A; Figure 5).

Stress Field Fluctuations Modulated by Seismic Cycling
Fluctuations in stress magnitudes can be caused by seismic cycling.It has been reported that earthquake events generate stress drops of 0.01-100 MPa, depending on the rheology, roughness of fault, geometry of slip area, and heterogeneous stress fields (Allmann & Shearer, 2009;Baltay et al., 2011;Candela et al., 2011;Cocco et al., 2016;Oth et al., 2010).The observation of localized normal faulting regimes in the HSM may be related to seismic cycling in the region.The normal faulting regimes observed along central HSM boreholes Kauhauroa-2 (1980-2075 m TVD), Kauhauroa-5 (1330-1345 m TVD), and Tuhara-1A (1700-1820 m TVD) occur where S Hmax and S v are very similar and are greater than S hmin (Figures 5, Figures 10b and 10c).In such stress state scenarios, a post-seismic stress drop of only a few MPa after great earthquakes or frequent moderate earthquakes in the HSM region could perturb the delicately balanced stress magnitudes surrounding these boreholes, switching σ 1 = S Hmax to σ 1 = S v that is, from a reverse/strike-slip to a normal stress state, accompanied by small rotations in the S Hmax orientation.Such rotations in borehole-derived S Hmax orientations and changes in the stress state from thrust faulting regime to a normal faulting regime have been reported following 1999 M w 7.6 Chi-Chi earthquake in Taiwan and 2011 M w 9.0 Tohoku-Oki earthquake in Japan (W.Lin et al., 2007Lin et al., , 2013)).

Study Limitations
Although the results of this study provide a comprehensive quantification of the three principal stress magnitudes within the HSM for the first time by integrating all available borehole data (from 44 boreholes), it is important to note that the results of S Hmax values from BO method are subject to limitations owing to implicit assumptions and a lack of certain data.
Estimations of S Hmax magnitudes can be significantly influenced by UCS values, especially when direct UCS data are unavailable and there are no established empirical relationships to estimate UCS from other rock properties in this region.Although we used a range of empirical relationships developed worldwide to include uncertainty associated with the calculated UCS values into the stress model, it should be noted that the lack of the regional data could result in an imprecise estimation of the S Hmax magnitude and therefore stress state within the HSM.
The magnitude and orientation of S Hmax are determined by the analysis of BO geometries at the time of drilling, assuming that BOs do not grow over time.However, studies reveal that in porous, fractured, and weak muddy intervals, new BOs can form or existing BOs can grow both along and radially around the borehole wall over time (Ask & Ask, 2007;J. C. Moore et al., 2011), thus introducing errors into inferences of stress orientations and stress magnitudes.
Furthermore, this method assumes that BOs occur when the circumferential stress exceeds the UCS.However, in reality, the borehole wall is only unconfined on one side.Therefore, the UCS represents the minimum threshold that the circumferential stress must exceed, and failure will actually occur somewhere between the UCS and the biaxial rock compressive stress.Also, it should be noted that this method assumes BOs form when circumferential stress near the borehole is the maximum stress; however there are features like knock-outs, which resemble BOs, are formed when the axial stress is the maximum stress (Bowes & Procter, 1997;Ding et al., 2022;Pašić et al., 2007).Given the depths of these wells, and HSM tectonic setting, it is likely that maximum circumferential stress exceeds axial stress, but it remains a potential issue with the BO width method.
In addition, there is uncertainty involved in measuring BO widths using resistivity image logs as these logs may not capture a full coverage of the borehole circumference, leading to incomplete imaging of BO widths (Kingdon et al., 2016;Massiot et al., 2019).Therefore identifying breakouts using such images may involve some degrees of subjectivity, which in turn can impact the estimation of S Hmax .
Another source of uncertainty in S Hmax calculation is the estimation of σ 3 from LOTs.In addition to the absence of LOT measurements along each borehole (Section 4.2), this study commonly uses LOP (representing the pressure required to initiate a fracture in the stress-disturbed near a borehole) obtained from LOTs to estimate σ 3 , as fracture closure pressure points and pressure-time plots are rarely available from XLOTs or hydraulic fracturing tests.However, studies suggest that LOP typically tends to slightly overestimate σ 3 (less than 2 MPa), likely leading to an underestimation of S Hmax magnitude (Addis et al., 1998;Bell, 1996;Breckels & van Eekelen, 1982;White et al., 2002).
In addition to the lack of direct P p measurement tests such as RFTs and MDTs, we use mud weight during drilling, equivalent circulating density (ECD), to calculate S Hmax .Not only ECD is slightly higher than real Pp but also considering that BOs can grow over time, it is possible that BOs form under lower static mud weight than ECD (Chang & Song, 2016).Therefore, using ECD as an estimate of annular pressure may lead to overestimation of S Hmax.
Density logs are typically not acquired within the shallow depth intervals of drilled boreholes.As a result, extrapolation methods are employed to determine a complete S v profile.However, these extrapolation methods introduce uncertainty, which can lead to uncertainties in determining S Hmax.
Although this study takes into account all these uncertainties and quantify them where possible, future studies benefit from more accurate and extensive data to build a reliable stress model for the HSM and to better constrain faulting regime.Despite these limitations, this study reveals a dominant thrust to strike-slip faulting regime throughout the entire HSM, consistent with the long-term contractional tectonic deformations resulting from the subduction of the Hikurangi Plateau beneath the North Island and the strike-slip faulting generated by the forearc rotation of the East Coast.

Conclusions
This work represents the first comprehensive determination of the in situ stress state of the HSM margin using available borehole data.We found a σ 3 :S v = 0.6-1 at depths above 650-700 m TVD, while σ 3 :S v = 0.92-1 below this depth interval along the HSM.Stress magnitudes calculated from borehole data indicate that the S Hmax :S v ratios ranging from 0.95 to 1.81 in the central HSM and 0.95-2.3 in the the southern HSM.These principal stress magnitudes result indicate a σ 1 = S Hmax and a thrust to strike-slip faulting regime across both the central and southern HSM.The pre-existing NNE/NE striking reverse faults along both the central and southern HSM infer that stress regime was initially in a contractional state such that σ 3 : S v = 1, σ 1 = S Hmax , and a dominant NW-SE S Hmax , consistent with NW-SE component Pacific-Australian plate motion.Taking the contemporary stress state of σ 1 = S Hmax and ENE-WSW S Hmax orientation and initial stress state into account in the central HSM, these observations suggest that the compressional regime has shifted from subparallel to perpendicular to the NW-SE Hikurangi convergence direction overtime in this region.Variation of the central HSM stress state over geological time may result from forces arising from Hikurangi forearc rotation either by itself or facilitated by the upper plate, shallow, high overpressures in the central HSM.Along-strike variation in slip behavior may also play a role by releasing stress overtime due to SSEs and frequent earthquakes, hence changing the stress state in the central HSM, while in the southern HSM, the modern WNW-ESE/NW-SE σ 1 (S Hmax ) remains subparallel to NW-SE Hikurangi convergence direction overtime, may reflect the interseismic locked nature of the plate interface.Overall, both lower overpressure and interseismic locking could contribute to making it difficult for southern HSM faults to undergo strike-slip in response to forearc rotation forces.Finally, stress determination highlights localized normal stress states within the HSM forearc interpreted to be due to processes such as the presence of localized active normal faults or fluctuations in stress magnitudes modulated by seismic cycles.The determination of HSM in situ stresses in this study will provide an invaluable information for improving our understanding of the stability of upper plate faults and will facilitate more quantitative efforts to assess the seismic hazard potential of the HSM that will support disaster risk reduction plans.This publication has emanated from research supported in part by grants from Science Foundation Ireland (SFI) under Grants 13/RC/2092 and 17/RC-PhD/3481 and is co-funded by Geological Survey Ireland (GSI).We appreciate the Editor Whitney Behr, Weiren Lin, and an anonymous reviewer for constructive and insightful reviews, which have greatly improved the manuscript.Special thanks to Dr. Mark Tingay for his comments on the limitations of this study.The authors would like to thank the New Zealand Petroleum and Minerals group (NZPM) within the Ministry for Business, Innovation and Employment (MBIE) for providing access to borehole data and supporting materials for this study.We also thank Schlumberger for providing academic licenses for Techlog 2018.1 to University College Dublin and University of Liverpool.We thank MathWorks for providing academic licenses for MATLAB to University College Dublin.We thank Esri for providing the academic license of ArcGIS Pro to University College Dublin.For the purpose of Open Access, the author has applied a CC BY public copyright license to any Author Accepted Manuscript version arising from this submission.

Figure 3 .
Figure 3. Map showing the distribution of the studied boreholes along the Hikurangi Subduction Margin.Red and pink circles represent boreholes with image logs and oriented four-arm caliper data, where S Hmax orientations determined by Behboudi et al. (2022).(b) Schematic cross-section across the Hikurangi subduction zone (line A-A′, (a)).The red lines show the location of boreholes with image log data relative to the subduction interface.Abbreviations: NIDFB, North Island Dextral Fault Belt; TVZ, Taupo Volcanic Zone.(b) modified from Barnes et al. (2010).

Figure 4 .
Figure 4. (a) Borehole schematic showing local principal stresses (σ θθ , σ zz , and σ rr ) at the borehole wall as a function of azimuth (θ) measured relative to S Hmax orientation and presence of breakouts for an example in which S Hmax = 50 MPa, S v = 45 MPa, S hmin = 40 MPa, and UCS = 45 MPa.The red shaded region schematically shows the circumference where σ θθ is large enough to exceed the compressional strength of the formation and induce BOs.(b) Diagram of a borehole cross-section showing the relationship between BOs, DITFs, and the horizontal principal stress orientations.(c) Example of DITFs as they appear on a resistivity image log.(d) Examples of BOs as they appear on a resistivity image log.Figures 2c and 2d are from resistivity image logs of borehole Kauhauroa-5.Abbreviations: UCS = unconfined/uniaxial compressive strength; BO = borehole breakout; DITF: drilling induced tensile fracture.

Figure 5 .
Figure 5. Calculated far-field in situ stress magnitudes with depth in borehole Tuhura-1A.(a) P-wave slowness (blue line) de-spiked and smoothed over 3 m intervals (pink line).(b) Range of UCS values derived from P-wave slowness using relations in Tables S4 and S5 in Supporting Information S1.(c) Azimuth of borehole BOs.(d) Calculated S v (solid black line), σ 3 (green field), and S Hmax (red field) magnitudes by considering that pore pressure is equal to mud weight.The hydrostatic pressure (gray line) is computed assuming a sea water density of 1.03 g/cc.The σ 3 and the range of σ 3 is determined from the average σ 3 :Sv = 0.95 and σ 3 :Sv = 0.92-1, respectively.(e) Tectonic stress regime index.Abbreviations: Δt C : P-wave slowness; UCS = unconfined/uniaxial compressive strength; BO = borehole breakout; FIT = formation integrity test; TF failure: thrust faulting failure; μ: friction coefficient; S v : vertical stress; S hmin : minimum horizontal stress; S Hmax : maximum horizontal stress; σ 3 : minimum principal stress; A ϕ : Tectonic stress regime index.

Figure 6 .
Figure 6.S v gradient profiles from boreholes in the (a) Hikurangi Subduction Margin (HSM) onshore and (b) HSM offshore.S v : vertical stress; TVD: true vertical depth from ground level for onshore boreholes and sea level for offshore boreholes.

Figure 7 .
Figure 7. Results of the leak-off test run at 722.5 m TVD in borehole Tawatawa-1.Pressure versus volume of mud pumped to the formation curve reveals that the leak-off pressure (LOP) is 12.7 MPa.

Figure 8 .
Figure 8. σ 3 : S v along the Hikurangi Subduction Margin.Abbreviations: LOT: leak of test; FIT = formation integrity test; S v : vertical stress; σ 3 : minimum principal stress; TVD: True vertical depth from ground level for onshore boreholes and sea level for offshore boreholes.

Figure 9 .
Figure 9. Graph showing Pp profiles for the central Hikurangi Subduction Margin (HSM) (blue lines) and the southern HSM (red lines) calculated from drilling mud weights.Direct pore pressure data are from repeat formation test (RFT) in borehole Kauhauroa-5 (central HSM) and modular dynamic test (MDT) in borehole Titihaoa-1 (southern HSM).

Figure 11 .
Figure 11.Analysis of stress magnitudes using stress polygon defined by Coulomb friction law with a friction coefficient (μ) of 0.6 in borehole Kauhauroa-5 at a depth of 1742.3 m where DITFs are observed.The green shaded area represents σ 3 values, ranging from the average σ 3 : S v = 0.95 to σ 3 : S v = 1.The blue shaded area represents S Hmax range which local hoop stress is tensile, and DITFs are formed.NF: normal faulting, SS: strike-slip faulting, TF: thrust faulting, DITF: drilling induced tensile fracture, μ: friction coefficient, S v : vertical stress, S Hmax : maximum horizontal stress, σ 3 : minimum principal stress.

Figure 13 .
Figure 13.Graph shows (a) p-wave velocity, (b) rock strength (UCS), and (c) Young's modulus in clay and sand-siltstone sediments as a function of depth across the central and southern HSM.

Table 1
Leak-Off Pressure (LOP) Measurements, Formation Integrity Test (FIT) Data, and S v for Boreholes in the Central and Southern HSM

Table 1 Continued
a True vertical depth from ground level for onshore boreholes and sea level for offshore boreholes.b Formation integrity test.c Leak-off pressure.d S v vertical stress.e LOP or FIT is equal or greater than S v .