Dynamical Modeling of Fault Slip Rates at the New Zealand Plate Boundary Indicates Fault Weakness

We construct a thin‐sheet dynamical model of the New Zealand plate boundary that includes faults. Our model fits fault slip rates, style of distributed deformation, and is constrained by relative plate boundary motion. We assume a pseudo‐plastic rheology and achieve a best fit to slip rate observations with a deviatoric stress magnitude of 20 MPa. Modeled local forces are significant at Puysegur and Hikurangi subduction zones, and smaller forces are related to mantle downwelling beneath South Island and Havre Trough mantle upwelling. Modeled tractions on faults are mostly 5–20 MPa, similar to or slightly smaller than stress magnitudes adjacent to faults. Modeled shear tractions are generally 2–10 MPa, comparable to stress drops during earthquakes. Modeled stress orientations and fault dips suggest that many faults are not optimally oriented for their style of faulting. Notably small traction magnitudes of <5 MPa and shear tractions of <0.5 MPa are modeled for faults in the central North Island Dextral Fault Belt (NIDFB), which we infer to be very poorly oriented. Friction coefficients on faults (ratio of shear stress to effective normal stress) are in the range 0.1–0.3 for major crustal faults such as the Alpine Fault and Marlborough faults, but subduction zones and the NIDFB have values <0.1. We propose that low values of long‐term fault strength, shear stress resolved onto the fault, and overall magnitudes of deviatoric stress in the crust are a consequence of dynamic weakening of faults during fault slip.


Introduction
Fault strength determines the magnitude of stress required to produce slip on the fault, for example, during an earthquake.With the strength of the upper crust often assumed to be limited by the frictional strength of faults (e. g., Brace & Kohlstedt, 1980;Bürgmann & Dresen, 2008), fault strength provides a constraint on stress magnitudes in the upper crust.However, it is difficult to measure fault strength directly and it is unresolved whether faults are generally weak or strong (e.g., Klein et al., 2009;Scholz, 1998;Townend & Zoback, 2000).We address this question through an analysis of regional faulting patterns in New Zealand.
Fault strength is often quantified in terms of the friction coefficient, which is defined as the ratio of shear stress to effective normal stress resolved onto the fault plane.Static friction coefficients for rocks measured directly in laboratory experiments are generally found to be 0.6-0.9(e.g., Byerlee, 1978;Kohlstedt et al., 1995;Rouet-Leduc et al., 2018).Fault strengths have also been estimated from stress magnitudes calculated in dynamical models of crustal deformation, which generally infer lower friction coefficients of <0.25 (e.g., Bird & Kong, 1994;Klein et al., 2009;Liu & Bird, 2002).Klein et al. (2009) indirectly estimated long-term average fault strengths in western North America by using modeled stress differences as a proxy within deforming regions.Liu and Bird (2002) incorporated faults into a thin shell model of New Zealand and found a uniform best fit fault friction coefficient of ∼0.17.Similar thin shell modeling has also found low friction coefficients of 0.1-0.25 for other plate boundaries (Bird & Kong, 1994;Bird et al., 2008).For the New Zealand plate boundary zone, deviatoric stress magnitudes have been estimated to be <35 MPa (Hirschberg & Sutherland, 2023;Hirschberg et al., 2018;Lamb, 2015), but these models treat the lithosphere as a continuum (in the context of a viscous thin sheet) and do not estimate fault properties directly.Hirschberg and Sutherland (2022) modeled faults in the New Zealand plate boundary zone, however, that model only investigated kinematics (including slip rates and distributed strain rates) and did not estimate stress, fault strength, or other dynamical properties.
We use the faulted thin-sheet method of Haines and Sutherland (2018) to model the dynamics of the New Zealand plate boundary, including fault strength and tractions on individual faults.We fit fault slip rates from the New Zealand Community Fault Model (NZCFM) v1.0 of Seebeck et al. (2023), Hikurangi subduction zone rates from a plate boundary kinematic model (Hirschberg & Sutherland, 2022), which also provides an estimate of the "missing" distributed deformation.Based on our preferred model of plate boundary rheology (Hirschberg & Sutherland, 2023), we assume a pseudoplastic, equal stress rheology with a best fit stress magnitude of ∼20 MPa, although we test sensitivity to the stress magnitude assumption.

Tectonics of the New Zealand Plate Boundary
At the New Zealand plate boundary, the Pacific and Australian plates are obliquely converging at ∼40 mm/yr (Altamimi et al., 2012;DeMets et al., 2010; see Figure 1).To the southwest of South Island (SI), the Australian Plate is subducting beneath the Pacific Plate at the Puysegur subduction zone.Through SI, ∼80% of plate motion is accommodated on the dextral-reverse Alpine Fault, which has a slip rate of ∼27 mm/yr (Howarth et al., 2018;Norris & Cooper, 2001).The remaining motion is accommodated in a wide band of faulting through southeastern and eastern SI (Seebeck et al., 2023).
At the northern end of the Alpine Fault, plate motion is partitioned onto the Marlborough Fault System (MFS).Hope Fault has the largest slip rate in MFS accommodating 23 ± 4 mm/yr (Langridge et al., 2003) and other major faults in MFS include Wairau, Awatere, and Clarence faults with slip rates of ∼4-6 mm/yr (Benson et al., 2001;Van Dissen & Nicol, 2009;Zachariasen et al., 2006).Reverse faulting occurs west of MFS in northwest SI (Seebeck et al., 2023).
East of North Island (NI), the Hikurangi subduction zone is the southern extension of the Tonga-Kermadec-Hikurangi subduction zone and involves the subduction of the thickened oceanic crust of the Hikurangi Plateau (Williams et al., 2013).In the overriding Australian Plate, dextral motion is accommodated in the North Island Dextral Fault Belt (NIDFB) (Beanland, 1995).In the southern NIDFB, the dextral-reverse faults include Wellington Fault with a slip rate of 6.8 ± 0.8 mm/yr (Van Dissen & Berryman, 1996) and Wairarapa Fault with a slip rate of 11.0 ± 0.5 mm/yr (Wang & Grapes, 2008).In the northern NIDFB, the sense of slip is predominantly dextral-normal (Mouslopoulou et al., 2007).

Method
We construct a dynamical thin-sheet model of the New Zealand plate boundary zone that includes faults and fits observations based on the method of Haines and Sutherland (2018).We fit observations of fault slip rates in a manner consistent with plate boundary motions and the style and direction of distributed deformation.We specify initial parameters representing strengths of faults, strengths of elements between faults, and forces acting within the model.We use these strength and force parameters to produce an "a priori" solution according to dynamical principles.We then calculate an "a posteriori" solution that combines the parameters with observations and  (Seebeck et al., 2023).Dashed blue contours indicate plate interface depth (in kilometers) for the Hikurangi subduction zone (Williams et al., 2013).MFS: Marlborough Fault System; NIDFB: North Island Dextral Fault Belt; TVZ: Taupō Volcanic Zone.
departs from dynamical principles in order to better fit observations.We use differences between the a priori and a posteriori solutions to adjust the parameters to improve the fit to observations.We iterate over the a priori solution, a posteriori solution, and parameter adjustment until an optimized fit to observations is achieved.

Dynamical Principles
We assume parameter values corresponding to material strength and local forces and use dynamical principles relating these parameters to calculate the a priori solution.Strain rate ėij in an element is related to stress τ ij by the locally linear relationship where the strain rate capacity η specifies the weakness of the material.Subscript indices refer to horizontal components only.
Similarly, slip rate [u i ] on a fault segment is related to horizontal traction on the fault τ jk n k by the relationship where the slip rate capacity κ ij specifies the weakness of the fault.A large slip rate capacity represents a weak fault that slips easily in response to an applied traction while a small slip rate capacity represents a strong fault that requires a larger traction to slip at the same rate.Haines and Sutherland (2018) show that κ ij can be expressed in terms of a single slip rate capacity κ and the geometry of the fault segment: where n i is the horizontal normal to the fault segment, θ is the dip of the fault, and δ ij is the Kronecker delta.n i is determined by the strike of the fault in the model whereas κ and θ can be varied with each iteration of the model.In general, the slip rate capacity of a fault is not known so we allow it to vary freely in the model (provided κ > 0).The dip of a fault through the crust is often poorly or only partially constrained so we allow θ to vary within bounds.Equation 2 can be expressed in terms of transverse and horizontal normal components of slip [u t ] and [u n ] as where κ c = κ cos 2 θ and t k is the unit vector in the transverse direction.
Forces f i originating within the model are included in the form of force potentials ϕ ij where ∂ϕ ij /∂x j = f i .Force potentials are related to stress as which is analogous to the force balance equations ∂τ ij /∂x j + f i = 0.
Equations 1-5 are solved to find the a priori solution by minimization of the work rate functional with where κ inv ik is the inverse of κ ik .E( ė) represents the contribution to the work rate from strain rate in elements, and F ([u]) represents the contribution from slip rate on fault segments.Equation 6 is minimized in terms of velocities at element vertices, with velocities specified on the boundaries of the model, holding ϕ, κ, and η constant.

Observation Fitting
We additionally solve for a separate a posteriori solution that departs from physical principles to better fit observations.This is achieved by minimization of the functional where ẑ is the model-predicted value, z is the observed value, and σ is the observation's standard error.Weighting between dynamical principles and observations is determined by the factor α. Small values of α produce solutions that better fit observations but depart more from dynamical principles compared to larger values of α.
We use a value of α where observations and dynamical principles are expected to have equal contribution to L POST .This is N model is the number of degrees of freedom in the model, equal to the number of unknown velocity quantities in the model.E 0 ( ė) and F 0 ([u]) are the values of E( ė) and F ([u]) from Equation 7with ϕ ij = 0 everywhere.

Parameter Adjustment
Differences between the a posteriori and a priori solutions indicate how force potentials ϕ ij and the material parameters of strain rate capacity η and slip rate capacity κ ij can be adjusted to better fit observations.The a posteriori solution provides a better fit to observations than the a priori solution.The principle is, therefore, to adjust parameters in order to move the a priori solution of the next iteration closer to the a posteriori solution of the previous iteration, thereby improving the fit to observations of the a priori solution in the next iteration.
To estimate the change in force potentials δϕ ij that will efficiently improve the fit to observations, we consider the change in strain rate in elements δ ėij resulting from δϕ ij or a change in strain rate capacity δη.The change δ ėij is where τ ij = τij + ϕ ij and τij represents the contribution to the stress from the velocity boundary constraints (see Section 2.5) on the model.Changes to τ are assumed to be small (δτ ≈ 0) and therefore neglected for the purposes of constructing a parameter adjustment method.By setting δ ė ≈ ėpost ėpre , where ėpre and ėpost are the strain rates in the a priori and a posteriori solutions respectively, the parameter adjustment will result in a new a priori solution that is close to the old a posteriori solution and therefore a better fit to observations than the old a priori solution.
When adjusting force potentials, we set δη = 0, simplifying Equation 10 to δ ėij ≈ ηδτ ij .The new force potential ϕ new in terms of the old force potential ϕ old is Journal of Geophysical Research: Solid Earth As κ c is related to κ by the dip of the fault, we constrain the ratio of κ c and κ so that the dip is within 20°of the initial value (and 0 < θ < 90°).

Pseudo-Plastic Rheology
There is a trade-off between material strengths (κ and η) and force potentials: doubling the strengths while halving the force potentials may result in the same solution.We constrain the solution by imposing a pseudo-plastic, equal stress rheology, which we have previously shown is appropriate for the New Zealand plate boundary (Hirschberg & Sutherland, 2023;Hirschberg et al., 2018).
To effect this constraint, we adjust η as where τ mag is the element's stress magnitude and τ tar is the target stress that the elements will be adjusted to have.Strain rate capacity in elements with τ mag > τ tar will be increased, making the element weaker and reducing the stress in the next iteration.The reverse is true for elements with τ mag < τ tar .As this adjustment is applied at the start of each iteration, stresses in the model solution will not exactly match τ tar but will in general be close to it.τ tar is uniform across the model and for all iterations in a model.We test several models with different values of τ tar (see Section 3.1).

Model Construction
We solve the model in terms of velocities on a finite element mesh with maximum mesh size of ∼30 km.Faults are included in the model as locations where velocity discontinuities occur and the velocity is solved for on each side of the fault at each point along the fault.We use the same fault geometry as HS22-CFM of Hirschberg and Sutherland (2022), which is simplified from the NZCFM v1.0 (Seebeck et al., 2023).
Velocities are fixed at boundaries of the model.The eastern and western boundaries are treated as rigid with velocities determined by the relative plate motion rates of the Australian and Pacific plates, respectively.We use MORVEL plate velocities (DeMets et al., 2010).Velocities on the northern and southern boundaries are determined by interpolating between the plate motions on either side.Extension of 20 mm/yr at an azimuth 135°is set for Havre Trough (Caratori Tontini et al., 2019;Wright, 1993).No deformation is applied elsewhere on the boundary except for Hikurangi subduction margin, which is assumed to accommodate the remainder of the relative plate motion.On the southern boundary, 14 mm/yr of dextral slip is set for Puysegur Ridge and the remainder of the relative plate motion is applied to the Puysegur subduction margin.

Observations in the Model
The primary observations in the model are slip rates from NZCFM (Seebeck et al., 2023).Parameters in the NZCFM, including slip rate, rake, and dip, were inferred from compiled field and marine observations, with uncertainties in the parameters assigned reflecting the quality and variation in the observations (see Seebeck et al., 2023 and references therein).NZCFM includes variations in these parameters for different sections along faults where data supports it.We convert the rake-parallel rates of NZCFM to horizontal slip rates, propagating the slip rate and dip uncertainties (Figure 2).NZCFM does not include slip rates on the Hikurangi subduction margin.Instead, we use slip rates and uncertainties from the HS22-CFM model of Hirschberg and Sutherland (2022).To reflect difficulty assigning deformation to different sections of the southern Hikurangi Subduction Zone, we assign the southern subduction margin an uncertainty equal to the shortening rate in this region.
The simplified mesh geometry does not include small or closely spaced faults.Instead, these faults are included in the model as strain rate observations.Slip rates are converted to strain rates using Kostrov summation.The strain rate is where l is the length of the fault segment and A the area of the element.We also include a uniform strain rate of 1.5 ± 0.4 × 10 7 /yr over the width of Havre Trough, equivalent to a total oblique extension rate of 20 ± 5 mm/yr (Caratori Tontini et al., 2019;Wright, 1993).
To ensure that modeled distributed deformation is consistent with the style of faulting observed on nearby faults, we include observations of strain rate style.We construct a model of strain rate style using the method of Hirschberg and Sutherland (2022) (Figure 2b).We spread fault slip rates across surrounding points using a 25 km triangular taper to create an equivalent strain rate according to Equation 16.Contributions from all faults are summed to create a strain rate field with values in all elements.We quantify the strain rate style in terms of A 0 , defined as the ratio of the dilatation component of strain rate to the strain rate magnitude (adapted from Klein et al. (2009)) , which represents isotropic horizontal expansion.A value of 1 represents uniaxial contraction (e.g., from parallel reverse faults), 1 represents uniaxial extension (e.g., from parallel normal faults), and 0 represents simple shear (e.g., from parallel strike slip faults).
Observations of strain rate style are included as a set of observations separate from the strain rate observations derived from small and/or closely spaced faults.We set the standard error of A 0 to ±1 for all elements to reflect that it is not our primary observations and that the strain rate style may not exactly match the style of faulting.The large uncertainty allows the strain rate style to deviate from the style of faulting if required by the dynamics.

Initial Parameters and Optimization
We use a uniform initial strain rate capacity of 3 × 10 24 Pa 1 s 1 and initial force potentials from the pseudoplastic PLAS model of Hirschberg and Sutherland (2023), with the areal component of these force potentials ranging from 50 to 20 MPa.We use initial slip rate capacities proportional to the slip rate observation for each fault segment.We use the preferred dips from NZCFM when calculating initial slip rate capacity components and allow fault dips to vary up to 20°from their NZCFM value.
We use an iterative adjustment strategy of alternating adjustments to produce an a priori solution that minimizes misfit to observations while adhering to an (approximately) pseudo-plastic rheology.On odd-numbered iterations, we adjust slip rate capacity based on Equation 14.On even-numbered iterations, we adjust force potentials according to Equation 11 and strain rate capacities according to Equation 15.We smooth the force potentials using a Gaussian filter with σ = 50 km.
The overall solution method is summarized as: 1. Setup: specification of observations of slip rates, strain rates, strain rate styles, and velocity boundary conditions; and specification of initial values for adjustable parameters of uniform η, κ proportional to slip rate observations, and ϕ from Hirschberg and Sutherland (2023).2. A priori solution: compute a solution based on dynamical physical principles (minimization of work rate), that is, Equation 6 is minimized in to find a velocity field while holding η, κ, and ϕ fixed.3. A posteriori solution: compute a posteriori solution that compromises between dynamical principles (minimizing work rate) and minimizing the misfit to observations, that is, Equation 8 is minimized to find a velocity field while holding η, κ, and ϕ fixed.4. Adjustment of κ: slip rate capacity is adjusted based on the difference between a priori and a posteriori solutions using Equation 14; and the ratio of slip rate capacity components κ c /κ = cos 2 θ is constrained so that the dip θ is within ±20°from its initial NZCFM value. 5. New solution: repeat steps 2-3 to calculate new a priori and a posteriori solutions.6. Adjustment of ϕ and η: adjustment of ϕ in elements based on differences between a priori and a posteriori solutions using Equation 11; and simultaneously adjust of η to approach a pseudo-plastic rheology using Equation 15. 7. Iteration: Iterate steps 2-6 until the a priori solution stabilizes on the optimal solution that minimizes misfit to observations while approaching a pseudo-plastic rheology.

Overview of Preferred Model
We test a range of target stress magnitudes for the pseudo-plastic model rheology between 5 and 50 MPa (Figure 3).We assess the models based on the residual sum of squares of the model's fit to the 412 fault slip rate observations (weighted by the observations' standard errors).We use fault slip rate observations to assess the model fit as they are our primary observation.There is a clear minimum in the residual sum of squares at ∼15-20 MPa, with the best fit to observations occurring for a stress magnitude of 20 MPa, which we use as our preferred model.This deviatoric stress magnitude is ultimately constrained by variations in topography and bathymetry from Hirschberg and Sutherland (2023) and is consistent with the <35 MPa previously estimated for the New Zealand plate boundary zone (Hirschberg & Sutherland, 2023;Hirschberg et al., 2018;Lamb, 2015).For comparison, we show the results for a target stress magnitude of 30 MPa, the preferred average stress magnitude of Hirschberg and Sutherland (2023), in Supporting Information S1.Our best-fitting 20 MPa dynamical model provides an acceptable fit to fault slip rate observations (Figure 4) with an average misfit of ∼1.4 per observation.The model also produces strain rate styles that generally match the model of strain rate style (Figure 5).One exception is the lack of a zone of almost simple shear in the region adjacent to Alpine Fault, where a larger contraction component is modeled for the region.
While strain rate directions were not fit to observations, modeled strain rate directions generally match those from Figure 2. The direction of maximum compressive strain rate is generally 100°-120°across SI.Strain rate directions are ∼135°in eastern NI east of NIDFB and 030°-050°west of NIDFB.As we assume an isotropic effective viscosity, the modeled strain rate direction is the same as the modeled stress direction.

Forces
Deformation in the model is driven by velocity boundary conditions, which represent far-field forces from outside the model domain, and gradients in force potentials, which represent forces originating within the model.Modeled force potentials include topographic force potentials, derived from variations in topography and bathymetry, and non-topographic force potentials, derived from other sources within the model domain such as basal tractions.We estimate nontopographic force potentials by subtracting the topographic force potentials from the total force potentials estimated by the model (Hirschberg & Sutherland, 2023).The general pattern of forces are consistent with tectonic expectations, which supports the assumption of a pseudo-plastic rheology (Hirschberg & Sutherland, 2023), and is confirmed by our new model that contains faults.
Forces originating within the model domain (Figure 6a) are largest adjacent to the Hikurangi and Puysegur subduction zones, which is expected from tractions on the subduction thrusts (Hirschberg & Sutherland, 2023).Forces (per unit volume) along eastern NI are 100-400 N m 3 in southeasterly direction, consistent with tractions on the interface of the Hikurangi subduction zone.Similarly, forces in southwestern SI and to the southwest of SI are 200-500 N m 3 in an easterly-to-southeasterly direction, consistent with convergence at the Puysegur subduction zone.
Forces associated with Havre Trough, north of NI, are generally <100 N m 3 and have a divergent pattern consistent with extensional faulting.Such a pattern of basal tractions is consistent with the expected backarc mantle upwelling and corner flow associated with Hikurangi-Kermadec subduction (Hirschberg & Sutherland, 2023).
Non-topographic force potentials are negative in central SI and eastern NI (Figure 6b), indicating mantle downwelling.Local forces are expected from tractions on the subduction thrusts, but our model also suggests that basal tractions of an opposite sense produce forces of order 50 N m 3 in eastern SI and central NI.The inference of mantle downwelling beneath SI is consistent  HIRSCHBERG AND SUTHERLAND with independent geophysical observations and inferred tectonic history (Pysklywec et al., 2002;Stern et al., 2000;Sutherland et al., 2000).

Effective Viscosity
Our model allows for the estimation of effective viscosity in elements, representing the bulk strength between the faults included explicitly in the model (Figure 6c).Because of the assumed pseudo-plastic rheology, regions of low effective viscosity correspond to regions of high strain rate.Low effective viscosities of ∼10 21 Pa s are modeled in Havre Trough, increasing southwards to ∼10 22 Pa s in the Taupo Volcanic Zone.Effective viscosities of 10 21 Pa s are also modeled in regions of geometrical complexity, such as fault bends or fault terminations, or other regions where high strain rates are required by strain compatibility.Effective viscosities of ∼10 22 Pa s are  consistent with our previous analysis and within the range expected, based on other analyses from other regions (e.g., England & Molnar, 2015;Flesch et al., 2001;Ghosh et al., 2013;Hirschberg & Sutherland, 2023).

Fault Strength
Fault strength is expressed in the model in terms of slip rate capacity κ, the ratio of slip rate to traction (Figure 7a); physically it represents the material weakness η integrated over the width of the fault zone.In general, the fastest slipping faults are modeled with the highest slip rate capacities.Slip rate capacities of ∼10 16 m s 1 Pa 1 are modeled for the major plate boundary faults: Alpine Fault, Puysegur subduction margin and northern Hikurangi subduction margin.Slip rate capacities of ∼10 17 m s 1 Pa 1 are modeled for other major faults including faults of the MFS and the southern Hikurangi subduction margin.Faults in Hikurangi and Fiordland margins generally have slip rate capacities of 10 18 -10 17 m s 1 Pa 1 and faults with small slip rates in northwestern and southeastern SI tend to have slip rate capacities of 10 19 -10 18 m s 1 Pa 1 .
Faults of the NIDFB are modeled with notably larger slip rate capacity compared to other faults with similar slip rates.This is especially true of the Ruahine and Mohaka faults in central NIDFB, where a slip rate capacity of ∼10 15 m s 1 Pa 1 is modeled, approximately an order of magnitude larger than other major plate boundary faults.The Ruahine and Mohaka faults are still modeled as weak when using a higher target stress of 30 MPa, but the smaller modeled slip rate capacities of ∼3 × 10 17 m s 1 Pa 1 are more comparable with other major plate boundary faults, that is, models with higher stress require them to have less anomalous weakness.

Fault Dip
Fault dip θ is inferred indirectly from the relative values of the two components of slip rate capacity with κ c = κ cos 2 θ for an isotropic slip rate capacity (Figure 7c).Inferred fault dips are nearly vertical (>80°) for most faults in MFS and NIDFB, which is mostly but not always steeper than the preferred dip of NZCFM that was used for the starting slip rate capacity.The inferred dip on central Alpine Fault is ∼60°, ∼10°steeper than NZCFM's preferred dip, while the inferred dip of ∼80°for southern Alpine Fault is barely changed from the preferred dip.
Inferred dips in the Hikurangi accretionary margin are shallow with a typical dip of ∼20°, generally shallower than NZCFM's preferred dips.The inferred dip on the Hikurangi subduction interface is very shallow, varying from ∼0°-20°.Inferred dips on the Puysegur subduction interface and in the Fiordland margin are similarly shallow at ∼0°-20°, shallower than NZCFM's preferred dips.
We split the fault segments into reverse, normal, and strike-slip faulting styles based on their modeled relative transverse and normal horizontal components of slip (Figure 8a).We exclude fault-terminating segments from this analysis as these are the segments most likely to be modeled with a different faulting style to their observed style.Of the remaining, non-terminating segments, the modeled style agrees with the observed style for 92% of segments.Strike-slip segments are predominantly near vertical with ∼50% of strike-slip segments modeled with an equivalent dip >80°, including ∼75% of strike-slip segments in NI.Reverse faults in NI predominantly have a shallow dip with ∼50% having an equivalent dip of 20°-30°, including most of the faults of the Hikurangi accretionary margin.Reverse faults in SI show a more even spread of dips.Normal faults, which are significantly less numerous in the model, are generally modeled with an intermediate dip with ∼50% having an equivalent dip of 40°-70°.

Fault Tractions
Modeled average tractions on faults are generally 5-20 MPa (Figure 6d).Tractions of 10-20 MPa are modeled for Alpine Fault and MFS.Smaller tractions of ∼10 MPa are modeled for Puysegur and Hikurangi subduction interfaces.Larger tractions of ∼20 MPa are modeled for faults in southeastern SI and faults in southern Hikurangi accretionary margin.Notably small tractions <5 MPa are modeled for the Ruahine and Mohaka faults, which also have anomalously large slip rate capacities.The remaining faults in the model generally have small tractions of 5-10 MPa.
When resolved onto the fault plane using the equivalent dip, modeled shear tractions are generally 2-10 MPa (Figures 7b and 8c).Shear tractions are ∼10 MPa for Alpine Fault and MFS and 2-7 MPa for Hikurangi subduction zone.Notably small shear tractions of <0.5 MPa are modeled for NIDFB.Temporarily neglecting body forces, modeled normal tractions resolved onto the fault plane are <20 MPa in both islands (Figure 8b).For both components, modeled tractions are generally smallest on normal faults at <5 MPa.Reverse faults in NI generally have small normal tractions of <4 MPa and transverse tractions of <8 MPa, whereas both components are generally larger in SI reverse faults.Transverse tractions on SI strike-slip faults are generally ∼3× larger than transverse tractions on NI faults.Normal tractions are similar for strike-slip faults in both islands.

Fault Orientation
The optimal orientation of fault planes can be predicted by calculating the plane where the shear traction most exceeds the resistance provided by the effective normal stress and fault friction (Anderson, 1905).For rocks with typical static friction coefficients of 0.6-0.9(Byerlee, 1978), this is predicted to be 25°-30°from the axis of maximum compressional stress.By assuming that one of the principal stresses is vertical, it is possible to predict the optimal strike and dip for different styles of faulting (Anderson, 1905).The optimal dip is predicted to be 60°-65°for normal faults, 25°-30°for reverse faults, and vertical for strike-slip faults.In the horizontal plane, the optimal strike is predicted to be 25°-30°from the orientation of the maximum horizontal compressive stress S Hmax for strike-slip faults, parallel to S Hmax for normal faulting, and perpendicular to S Hmax for reverse faulting.The angle between the optimal fault plane and S Hmax is dependent on the friction of the fault plane.A higher friction coefficient of 1 predicts this angle to be 22.5°and no friction predicts an angle of 45° (Anderson, 1905).Additionally, it is possible for faults to slip even when they do not follow Anderson's theory.This can occur if the effective normal stress is reduced by pore fluid pressure or if the principal stresses are not vertical and horizontal (Sibson, 1985).Both situations are likely true in the forearcs of subduction zones.
We estimate how well faults are oriented in the vertical plane using our inferred equivalent dips (Figure 8a) and modeled faulting styles.Almost all strike-slip faults are modeled with dips near the optimal dip of 90°.Reverse faults in NI generally dip 20°-40°, near the optimal dip of 25°-30°, but many reverse faults in SI are steeper, indicating that they are less well oriented, likely due to some being reactivated normal faults.Normal faults are generally modeled with equivalent dips of 50°-70°, near the optimal dip of 60°-65°.
We compare the strike of fault segments to modeled S Hmax averaged from the two elements either side of the segment (Figures 8d and 9).The strikes of the Alpine Fault and faults of the MFS are generally 30°-60°from S Hmax , indicating that only some of these faults are well oriented for strike slip.Strikes of reverse faults in southeastern SI, northwestern SI and the Hikurangi accretionary margin are generally >70°from S Hmax .Northern and central Hikurangi interface are also modeled with strikes >70°from S Hmax , whereas the southern Hikurangi interface, where a larger dextral component is predicted (Hirschberg & Sutherland, 2022), is modeled with strikes ∼60°from S Hmax .
Many faults of the NIDFB are modeled as being poorly oriented in the horizontal plane.Faults in southern NIDFB generally strike approximately parallel to S Hmax , while those in northern NIDFB generally strike perpendicular to S Hmax .While this is well-oriented for their secondary senses of reverse and normal slip, respectively, it is poorly oriented for their dominant strike-slip sense (and the fault dip is poorly oriented for reverse and strike-slip movement).Consequently, these faults have low shear tractions predicted and are required to be anomalously weak.This is comparable to the San Andreas Fault where high angles between the fault and the regional stress field have been observed (Lundstern & Zoback, 2020;Townend & Zoback, 2004;Yang & Hauksson, 2013) and laboratory measurements of fault gouge indicate that the fault is weak in some places (Lockner et al., 2011).
Although faults of NIDFB may be poorly oriented in relation to the regional stress field, they may be more optimally oriented in relation to local stress fields.The element size and force potential smoothing in our model means the stress field resolution is likely on the order of several 10 s of km and variation close to faults on length scales smaller than this cannot be resolved; and might be expected, given the local changes modeled near there.
The NIDFB is modeled as the boundary between the region east of the faults with S Hmax oriented ∼135°( associated with extension in the TVZ) and the region west of the faults with S Hmax 030°-050°(associated with compression along Hikurangi subduction zone), indicating that the stress field undergoes significant rotation across the NIDFB.

Fault and Lithosphere Strength
Vertical stress and differential stress in the uppermost portion of the crust increases approximately linearly with depth (e.g., Brudy et al., 1997), but direct measurement of stress is mostly limited to depths <5 km and there is uncertainty at greater depths.Lithospheric strength is thought to be dominated by the brittle strength of rocks at depths 10-20 km and by ductile strength that follows power-law flow at greater depths (Brace & Kohlstedt, 1980;Kohlstedt et al., 1995).The presence of weak zones, that is, faults, will alter the stress at which failure can occur and therefore the overall strength of the lithosphere.
Vertical stress at a given depth is simply the weight (per unit area) of the rock above it.In reverse faulting regimes, the minimum principal stress is vertical and therefore the mean stress is larger than the vertical stress.In normal faulting regimes, the maximum principal stress is vertical and therefore the mean stress is smaller than the vertical stress.This implies that the mean stress in reverse faulting regimes is higher than the mean stress in normal faulting regimes and hence that differential stress magnitudes at equal depth are expected to be dependent on the stress regime (Sibson, 1974).
The effective normal stress on a fault plane is reduced if pore fluid pressure is increased, reducing the shear stress required for failure (Hubbert & Rubey, 1959).Pore fluid pressure may vary transiently during the earthquake cycle, altering the effective normal stress and therefore the strength of the fault during the earthquake cycle (Sibson, 1994).
Fault strength may be expressed, according to Amontons's law (Amontons, 1699) and neglecting cohesion, in terms of the friction coefficient μ, defined as the ratio of shear stress magnitude |σ s | to effective normal stress σ n resolved onto the fault.In terms of modeled deviatoric stress τ and assuming that one of the principal stresses is vertical, the magnitudes of shear and effective normal stresses on the fault plane are where θ is the fault dip, λ is the ratio of pore pressure to lithostatic stress, and σ zz is the vertical stress.Assuming a constant density, σ zz = ρgz.Averaging over a layer of thickness h, the vertically averaged friction coefficient μ is As we use thin-sheet modeling, the modeled values are vertically averaged over the thickness of the sheet.We do not explicitly set the thickness of the sheet, but our fitting of surface observations assumes those observations are representative over the thickness of the sheet with no decoupling between the surface and the base of the sheet.Our introduction of local force potentials, which include the effects of basal tractions, means that we can model coupling between the base of our sheet and material beneath it.A pseudo-plastic rheology supports surface observations being related mainly to the strength of the crust (Hirschberg & Sutherland, 2023) and we interpret our results on the basis of our thin sheet representing the crust but we discuss below the implications of assuming different thicknesses.
Because we model in two dimensions, faults in the model are reduced to one dimension.Fault dip is represented as part of their anisotropic slip rate capacity and we assume no differences across the vertical profile of the fault.We do not account for changes in fault geometry with depth such as intersections of parallel faults or the possibility of faults widening at depth into distributed shear zones.
Notably small μ< 0.01 is modeled for Puysegur subduction zone and Hikurangi subduction zone is modeled with variable μ ranging from 0.01 to ∼0.1.Smaller μ are modeled near bends in the Hikurangi fault trace; this effect is related to our rheology assumption which does not allow stress to accumulate in zones of geometrical complexity, causing the model to compensate by locally weakening the fault and surrounding elements.
The Alpine Fault and MFS are modeled with μ of 0.1-0.3.Values of μ of 0.05-0.2are modeled for most other faults, including most in southeastern SI, northwestern SI, and in the Hikurangi accretionary margin.
Notably small μ < 0.01 is modeled for many faults of the central NIDFB.Larger μ of 0.05-0.1 is modeled at the northern end where it meets the TVZ and μ of 0.03-0.07 is modeled for Wairarapa Fault in the southern NIDFB.
The model with a higher target stress of 30 MPa estimates larger (but still low) μ of 0.01-0.03for the central NIDFB.Low estimates of μ are a consequence of the NIDFB being poorly oriented relative to the modeled stress (see Section 4.1), resulting in small shear tractions resolved onto the faults.This requires the faults to be weak in order to fit observed slip rates.
With the exception of anomalously low friction coefficients in the NIDFB, there is no clear difference in friction coefficients between normal, reverse, or strike-slip faults (Figure 10b).Friction coefficients on strike-slip faults in SI and reverse faults in both islands are approximately evenly distributed within the range 0-0.2.

Thin Sheet: Lithosphere or Brittle Crust
Our low values of μ are consistent with values inferred from previous dynamical models that solve for long-term averages.Liu and Bird (2002) modeled the New Zealand plate boundary using a best fit uniform effective friction of ∼0.17.Similar μ of <0.25 has been modeled for other plate boundaries (e.g., Bird & Kong, 1994; Bird  et al., 2008;Klein et al., 2009), suggesting that broader conclusions about the mechanics of the lithosphere and crust may be drawn from our results.
We use a layer thickness of 20 km when estimating μ, but other layer thicknesses also predict small friction coefficients.Using h = 100 km, representing a lithospheric thickness, decreases our inference of μ by a factor of ∼5, reducing the median μ from ∼0.1 to ∼0.02.This would indicate that fault friction is very small and much smaller than estimated by previous dynamical modeling (e.g., Bird & Kong, 1994;Klein et al., 2009;Liu & Bird, 2002).Using a thinner layer of h = 10 km increases μ by a factor of ∼2, increasing the median μ to ∼0.2.
Our modeled forces are consistent with the base of the seismogenic crust being the base of our thin sheet.Modeled non-topographic forces are interpreted as basal tractions of up to 15 MPa at the base of a 20 km thick layer.The basal tractions are comparable in magnitude to the modeled deviatoric stress and tractions on faults.The directions of modeled forces are consistent with those expected from lower lithosphere movements of convergence and subduction and/or downwelling acting on the base of the seismogenic crust, supporting the choice of the seismogenic crust as an appropriate sheet thickness.
For a 100-km thick layer (sheet), the forces implied are equivalent to much larger basal tractions of up to 75 MPa, which is inconsistent with our inference of stress magnitudes within the sheet and is inconsistent with the hypothesis that the base of the lithosphere is weak.Hence, we suggest that the strength peak of the brittle-ductile transition may act as the main strength guide within the lithosphere for determining surface deformation.We suggest that thin-sheet models based on power-law flow potentially overstate the strength of the upper mantle, possibly due to the presence of weak shear zones (e.g., Bürgmann & Dresen, 2008;Norris & Toy, 2014).
Finally, a note of caution about the limitations of the thin-sheet approach.The two-dimensional geometry of our model cannot fully capture the three-dimensional geometry of the New Zealand plate boundary.Many faults in the Hikurangi accretionary margin and the NIDFB likely terminate against the Hikurangi subduction interface (Henrys et al., 2013;Seebeck et al., 2023;Williams et al., 2013).This limits the depth over which vertical averages on these faults can be meaningfully considered to depths shallower than the subduction interface (which supports our choice of layer thickness).Model resolution may be an issue in the vicinity of the NIDFB because it occupies a boundary between different strain rate styles (Figure 5d) and hence the directions of principal stresses change dramatically over short horizontal distances comparable to the scale of our model discretization.
Changing the pore fluid pressure to increase our modeled friction coefficients by a factor of 5 would require a very wet crust with λ ≈ 0.9.Friction coefficients measured for some fault gouge materials are lower than for most other rocks, however, these are still generally 0.2-0.5 (e.g., Boulton et al., 2014Boulton et al., , 2017;;Byerlee, 1978;Carpenter et al., 2011;Morrow et al., 2007), the lower end of which is consistent with only our strongest modeled faults.
Our estimates of vertically averaged friction coefficients represent long-term values averaged over many earthquake cycles.Slip is localized on faults because they are weaker than the surrounding crust and therefore fail first when stress accumulates, allowing surrounding crust to maintain greater strength (e.g., Townend & Zoback, 2000).
The dynamic strength of faults during slip is notably less than their static strength, and several different weakening mechanisms may be responsible (Di Toro et al., 2011).Noda et al. (2009) show theoretically, given reasonable assumptions of material properties, that the thermal pressurization mechanism of weakening (Rice, 2006;Sibson, 1973) leads to earthquake rupture propagation at ratios of shear traction to normal traction (μ) of 0.2-0.3 and earthquake stress drops of ∼3 MPa, even though static friction coefficients are 0.6-0.9around the rupture front.Noda et al. (2009) also show that the average shear traction doing work during an earthquake is <10-20 MPa.
Earthquake data in New Zealand are generally consistent with our inferences of shear tractions resolved onto faults.Subduction-related earthquake stress drops near Gisborne and Puysegur had low values <0.2 MPa, whereas the Canterbury earthquake sequence in eastern SI had stress drops of ∼5-15 MPa (Fry & Gerstenberger, 2011).
The 1855 Wairarapa earthquake, which occurred on the southern end of the NIDFB, is the largest historic earthquake in New Zealand and perhaps the most enigmatic.The surface strike-slip rupture of ∼16 m is one of the largest recorded globally, and the short (∼145 km) length of rupture suggests an unusually large static stress drop (Rodgers & Little, 2006).Geomorphic evidence suggests this pattern of rupture has been repeated 7 times, so it was not an isolated anomaly (Manighetti et al., 2020).Our model predicts an anomalously low shear traction resolved onto the Wairarapa Fault, but it is probably one our least reliable results because: (a) stress orientations and slip rate style change rapidly across that region; and (b) the Wairarapa Fault connects and splays off the subduction thrust at about 20 km depth beneath Wellington (Henrys et al., 2013), so it cannot be considered in isolation.Wairarapa Fault ruptures likely start as subduction thrust ruptures and then propagate onto the Wairarapa Fault (Darby & Beanland, 1992).Moderately large (∼10 MPa) basal tractions are inferred on the subduction thrust by our model in this region (Figure 6b), and this may provide a viable explanation for the large surface ruptures of these anomalous earthquakes.
In summary, we suggest that our results are consistent with a variety of observations and previous models.The apparently low values of long-term tractions and fault strengths that we compute are caused by: weak fault gouges; dynamic weakening during earthquake slip; and elevated pore fluid pressure in the vicinity of subduction zones.Earthquake stress drops in New Zealand and globally (Allmann & Shearer, 2009) are consistent with the range of fault shear tractions that we compute.Heat flow data from the San Andreas Fault and the lack of notable temperature anomalies after earthquakes also require relatively small (<10 MPa) shear tractions resisting fault motion (Fulton et al., 2013;Kano et al., 2006;Lachenbruch & Sass, 1980;Lachenbruch et al., 1995;Williams et al., 2004).

Conclusions
We present a model of fault strength and of bulk strength between major faults for the New Zealand plate boundary zone using a pseudo-plastic, equal-stress bulk rheology.Known topographic gradients result in stress magnitudes and directions that are computed, and plate boundary motion is known and in a different direction, so an optimized solution is found that balances tectonic versus topographic forces to fit the spatial pattern of fault slip rate observations.We fit fault slip rate directions and rates from NZCFM (Seebeck et al., 2023).Our best fitting model has a mean deviatoric stress magnitude of 20 MPa, consistent with previous estimates for New Zealand (Hirschberg et al., 2018;Lamb, 2015).We model effective viscosities of 10 21 -10 22 Pa s for regions of the plate boundary that are rapidly deforming through distributed deformation.Modeled basal tractions are small except for those associated with subduction at Puysegur and Hikurangi subduction margins, mantle downwelling in central SI and eastern NI, and upwelling beneath Havre Trough.
Modeled horizontal tractions on faults are generally 5-20 MPa, comparable to the stress magnitude of ∼20 MPa modeled for surrounding material.Modeled shear tractions on faults are generally 2-10 MPa, comparable to stress drops in earthquakes (e.g., Allmann & Shearer, 2009;Boyd et al., 2017;Oth & Kaiser, 2014).Modeled stress orientations and fault dips suggest that many faults are not optimally oriented for their style of faulting.Notably small shear tractions of <0.5 MPa and correspondingly large slip rate capacities are modeled for faults in the central NIDFB, indicating that these faults are very weak.These faults are modeled with a regional stress field that is very poorly oriented for fault slip, but the regional stress field changes orientation significantly across these faults.
We estimate vertically averaged fault friction coefficients from our modeled fault tractions to be <0.2.These are much smaller than those measured from laboratory measurements of static rock friction (Byerlee, 1978), but consistent with low coefficients of friction previously estimated from plate boundary dynamical models (e.g., Bird & Kong, 1994;Klein et al., 2009;Liu & Bird, 2002).This suggests that the strength of faults during slip is notably less than the static strength of undamaged rock surrounding them.Our fault friction coefficients are estimated from fitting fault slip rates and represent long-term effective fault strengths over many earthquake cycles.Slip occurs on faults because they are weaker, and our estimates of strength are a type of average of resistance during slip, rather than static strength.The presence of faults reduces the strength of the crust, and limits crustal stresses in actively faulting regions to stresses that are much smaller than might be expected from Byerlee's law extrapolated to depths of 10-20 km (e.g., Brace & Kohlstedt, 1980;Bürgmann & Dresen, 2008).

Figure 1 .
Figure 1.Tectonic setting of the New Zealand plate boundary zone.Arrows indicate motion of Pacific Plate relative to Australian Plate (DeMets et al., 2010).Thick red lines indicate faults included in the model directly as localized slip rate and thin brown lines indicate faults included in the model indirectly as distributed strain rate; faults are from the New Zealand Community Fault Model v1.0(Seebeck et al., 2023).Dashed blue contours indicate plate interface depth (in kilometers) for the Hikurangi subduction zone(Williams et al., 2013).MFS: Marlborough Fault System; NIDFB: North Island Dextral Fault Belt; TVZ: Taupō Volcanic Zone.

Figure 2 .
Figure 2. Observations used in the model.(a) Slip rate observations on large faults are shown as colored lines and derived from New Zealand Community Fault Model(Seebeck et al., 2023), except for Hikurangi subduction margin, which is fromHirschberg and Sutherland (2022).Strain rate observation magnitudes derived from small faults (with the addition of Havre Trough) are indicated by blue-shaded polygons.(b) Model of strain rate style derived from the style of faulting on nearby faults and defined by Equation17.Yellow represents extension, blue is contraction, and white is simple shear (strike-slip faulting).Modeled directions of maximum compressional horizontal strain rate are indicated by black bars, however, these directions are not fit as observations during optimization of the model.

Figure 3 .
Figure 3. Residual sum of squares of the model's fit to the 412 fault observations (weighted by the observations' standard errors) for target average stress magnitudes of 5-50 MPa.The best fit to observations occurs for a target stress of 20 MPa with an average misfit of ∼1.4 per observation.

Figure 4 .
Figure 4. Misfit to fault slip rate observations for the preferred model.Small white circles indicate the location of the observation in the model.The misfits to observations (relative to each observation's standard error σ) are proportional to the diameters of the light blue circles.

Figure 5 .
Figure 5. Modeled kinematics of preferred model.(a) Modeled kinematic velocities relative to Australian plate.(b) Modeled kinematic velocities relative to Pacific plate.(c) Modeled strain rate magnitude.Largest strain rate magnitudes are modeled of 1-2 × 10 7 /yr for Havre Trough.Large strain rates are also modeled in regions of geometrical complexity as required by strain compatibility, including the southern Hikurangi-MFS transition.(d) Modeled strain rate style (colors) and directions (bars).Extension (yellow) is modeled for northern and western North Island and compression (blue) or oblique compression is modeled elsewhere.

Figure 6 .
Figure 6.(a) Modeled forces and force potentials from all sources within the model domain.The areal component of force potentials is indicated in color.Arrows indicate forces per unit volume.The largest forces are modeled adjacent to the Puysegur and Hikurangi subduction zones.Regions with negligible deformation where forces are not considered reliable are indicated in gray.(b) Modeled forces and force potentials not derived from variations in topography and bathymetry.(c) Modeled effective viscosity in elements (i.e., excluding faults).(d) Modeled stress magnitudes (colors) and modeled tractions in the horizontal plane (arrows).Stress magnitudes and tractions are generally near the target stress of 20 MPa.

Figure 7 .
Figure 7. Modeled slip rate capacities, shear tractions, and equivalent inferred dips on fault segments.(a) Slip rate capacity magnitude κ.Largest slip rate capacities are generally modeled for the fastest slipping faults.Notably large slip rate capacities are modeled for faults in the North Island Dextral Fault Belt.(b) Modeled shear tractions resolved onto the fault.(c) Equivalent dip implied by the relative values of the slip rate capacity components κ and κ c .(d) Difference between equivalent dip from the modeled slip rate capacity and preferred dip in New Zealand Community Fault Model (NZCFM)(Seebeck et al., 2023).The modeled equivalent dip was allowed to vary up to 20°from NZCFM's preferred dip.

Figure 8 .
Figure 8. Cumulative plots of equivalent dips and tractions for all non-terminating fault segments (i.e., excluding fault-terminating segments) and separately for reverse, strike-slip, and normal fault segments.Orange dashed lines indicate the proportion of fault segments, weighted by length, in North Island with that property below the specified value, green dotted lines indicate fault segments in South Island, and solid blue lines indicate all fault segments in the model.Normal fault segments are not separated by island due to the small number of normal fault segments.(a) Modeled equivalent dips.(b) Modeled magnitude of normal component of traction resolved onto the fault plane using the modeled equivalent dip and neglecting body forces.(c) Modeled magnitude of shear traction resolved onto the fault plane.(d) Difference between fault strike and modeled S Hmax .

Figure 9 .
Figure 9. (a) Difference between fault strike and orientation of modeled direction of maximum horizontal stress S Hmax .As we use an isotropic viscosity, S Hmax directions are the same as modeled directions of maximum horizontal strain rate shown in Figure 5d.(b) Difference between fault strike and the optimal fault orientation relative to S Hmax for each segment's modeled style of faulting.

Figure 10 .
Figure 10.(a) Modeled long-term vertically averaged fault friction coefficients μ estimated using Equation 19.Alpine Fault and Marlborough Fault System are modeled with μ of 0.1-0.3.Hikurangi is modeled with variable μ of 0.01-0.1.Notably small μ of <0.01 is modeled for faults of the North Island Dextral Fault Belt.(b) Cumulative plots of μ for non-terminating fault segments weighted by segment length and separated by fault style as per Figure 8. Orange dashed lines indicates North Island segments, green dotted lines South Island segments, and solid blue lines all segments.