A Kinematic Model of Quaternary Fault Slip Rates and Distributed Deformation at the New Zealand Plate Boundary

We construct a kinematic model of Quaternary long‐term deformation at the New Zealand plate boundary, including slip rates on major faults and strain rates between faults. We use an iterative method based on statistical and physical principles to find a velocity field that fits fault slip rate observations, has consistent off‐fault strain‐rate style, and is constrained by known plate motion. The kinematic model balances on‐fault and off‐fault deformation and provides improved estimates of fault slip rates and uncertainties in a statistically rigorous manner. We predict shortening rates of 45 ± 8, 34 ± 3, and 20 ± 3 mm/yr across the northern, central, and southern portions of the Hikurangi subduction margin, respectively; these rates are lower and better constrained by regional kinematics than previous estimates. In addition, the model predicts large strains adjacent to the Alpine Fault, indicating ∼9 mm/yr of plate motion in central South Island on faults with currently unknown slip rates. Differences between our long‐term velocity field and a contemporary velocity field arise mainly through interseismic locking on major faults. However, we suggest that differences with contemporary deformation in northeastern North Island are due to uncertainty in Havre Trough extension parameters and <5–10 mm/yr of missing dextral slip on unknown onshore faults, which implies seismic hazard is under‐predicted in that region.

Other kinematic models treat plate boundaries as a collection of fault-bounded blocks that can rotate but do not significantly deform (e.g., Avouac & Tapponnier, 1993;Peltzer & Saucier, 1996).While these rotating block models include faults, they do not account for deformation occurring between faults and, therefore, underestimate deformation.Block models may include geodetic velocities in their inversions, either exclusively or in combination with geological fault observations (e.g., McClusky et al., 2001;Wallace et al., 2004Wallace et al., , 2012;;Wang et al., 2017).However, the difference between the ∼10-year-timescale of geodetic observations documenting strain accumulated over only a portion of the seismic cycle and the >10,000-year-timescale of geological observations documenting strain release averaged over multiple earthquake cycles means that geodetic and geological rates can differ (e.g., Friedrich et al., 2003).Therefore, geodetically constrained kinematic models may better serve as a comparison to geologically constrained kinematic models than as a constraint on them.
Kinematic models are an essential input into many dynamic models investigating the forces and material properties observed in the kinematic model (e.g., England et al., 2016;Houseman & England, 1986;Medvedev & Podladchikov, 1999;Wdowinski et al., 1989).Thin-sheet models investigate vertically averaged viscosity and levels of stress in the crust or lithosphere, with many inverting for the dynamical model from the kinematic model, but they generally do not include faults in the model (e.g., Flesch et al., 2001Flesch et al., , 2007;;Ghosh et al., 2009;Hirschberg et al., 2018;Lamb, 2015).While not including faults in their thin sheet model, Klein et al. (2009) estimated fault strength by assuming that, long-term and depth-integrated, it is equal to stress differences in the seismogenic crust.Bird and Piper (1980) include the San Andreas fault as a narrow, weak zone in their forward model, but this was the only fault included in the model and they were only able to compare their model predictions to a limited set of observations available at the time.
We present a kinematic model of New Zealand that fits observations of deformation within the plate-boundary zone and is constrained by known plate-boundary motion, and our method could be applied elsewhere.We adopt the thin-sheet solution method of Haines and Sutherland (2018), which accounts for discontinuities at faults and does not require arbitrary grouping of observations.This non-linear regression method with plate-motion constraint uses physical principles to interpolate between observations.It is outside the scope of this paper to infer anything about absolute stresses in the crust or the relative strengths of faults, but we produce a fault slip-rate and crustal strain-rate field that is physically plausible, that is, off-fault strain rate could be generated by the same Through central South Island, the Alpine Fault accommodates 70%-80% of relative plate motion (Barth et al., 2014;Norris & Cooper, 2001).North of the junction with the Hope Fault, the Alpine Fault has a strike-slip rate of 13.6 ± 1.8 mm/yr and a reverse dip-slip rate of 3.4 ± 0.6 mm/yr (Langridge et al., 2010).South of this junction, the strike-slip rate increases to ∼27 mm/yr and the dip-slip rate to ∼7 mm/yr (Norris & Cooper, 2001).At the southern end of the Alpine Fault, from Fiordland to where it connects to Puysegur Subduction Zone, the Alpine Fault has a strike-slip rate of ∼23-31 mm/yr (Barnes, 2009;Barnes et al., 2005;Sutherland et al., 2006).Shortening in this region is mostly accommodated in the Fiordland accretionary wedge which shortens at 1-5 mm/yr (Barnes et al., 2002).
Southeast of the Alpine Fault is a broad zone of shortening and active mountain building.This zone is ∼100 km wide in Canterbury (central South Island) and ∼200 km wide in Otago to the south, with the difference likely due to variations in underlying lithosphere rheology (Upton et al., 2009).The total shortening accommodated on faults across this zone is ∼4 mm/yr, although this is likely to be a minimum in Canterbury due to the presence of faults with unknown slip rates (Litchfield et al., 2014).
To the south of South Island, subduction of the Australian Plate under the Pacific plate occurs at the Puysegur Subduction Zone.Earthquake slip vectors, seismic reflection images and seabed mapping indicate that the subduction thrust accommodates most of the plate motion (Melhuish et al., 1999).Puysegur Ridge Fault lies above the subduction thrust and is inferred to be strike slip but has an unknown slip rate (Collot et al., 1995;Delteil et al., 1996).
Early kinematic models of the New Zealand plate-boundary zone utilized repeated triangulations to calculate shear-strain rates and thus velocities on decadal timescales (Reilly, 1990;Walcott, 1978Walcott, , 1984)).The application of Global Positioning System (GPS) measurements has provided detailed knowledge of contemporary crustal deformation for onshore New Zealand (e.g., Beavan & Haines, 2001;Beavan et al., 2016;Haines & Wallace, 2020).GPS measurements only exist onshore but have been combined with geological data to create elastic block models that extend offshore (Wallace et al., 2004(Wallace et al., , 2007(Wallace et al., , 2012)).These models assume that contemporary and geological deformation are comparable, and conclusions are limited by the assumed geometry of blocks.Lamb (2015) applied the method of Lamb (2000) to fault slip rate and paleomagnetic data to invert for long-term velocities at 32 vertices.This kinematic model shows the first-order deformation pattern, including extension in the TVZ, shortening in the Hikurangi subduction zone and shear along the Alpine Fault and MFS.However, Lamb's model is calculated from a sparse mesh that treats many faults collectively rather than separately.His model, therefore, does not show detailed deformation features.

Model Objectives and Assumptions
We construct a kinematic deformation model across the New Zealand plate boundary zone that fits observations.This includes fitting slip rates on major faults and modeling distributed deformation between faults.
We apply relative velocities of the Pacific and Australian plates to the eastern and western boundaries of the deforming zone and treat these velocities as an exact constraint.We interpolate velocities between the two plates Sedimentary basins surrounding New Zealand show evidence of negligible deformation under the current tectonic environment and hence we identify locations of the eastern and western edges of the deforming zone based on previous mapping of seismic reflection data (Cook et al., 1999;King & Thrasher, 1996;Nathan et al., 1986;Sutherland et al., 2009).We treat the following regions as rigid: the undeforming Pacific plate to the east of New Zealand, the undeforming Australian plate to the west of New Zealand, and the Raukumara Basin northeast of North Island (see Figure S1 in Supporting Information S1). is the horizontal normal component of slip, with extension positive and shortening negative.Pred. is the value predicted by our best-fitting model.All slip rates in mm/yr.Slip rates for all faults are given in the supplementary information.We fit fault slip-rate observations (Figure 2), weighted by their uncertainties so that priority is given to fitting well-constrained observations (a maximum likelihood method).Most fault locations are put into the model domain (mesh; Figure S1 in Supporting Information S1) as a continuous series of discrete element edges along which velocity discontinuities are allowed.Small closely spaced faults are not entered as discrete mesh elements with velocity discontinuities (faults).Instead, observations of slip on each small fault are summed and converted into composite observations of the local strain rate within the triangular mesh element they occupy.Multiple slip or strain rate observations are allowed for the same mesh element.
Our model expects distributed deformation between faults, but we wish to obtain a physically meaningful model, that is, the on-fault and off-fault deformation should be of similar style (i.e., normal, reverse, strike-slip, or oblique) as nearby faults.To illustrate this necessity, consider a convergent plate boundary where the estimated fault slip rates sum to more than the plate motion.If the off-fault deformation style were free to vary in any way, the model would reconcile the excess summed convergent slip rates by distributing extensional strains between the faults, which is not physically plausible.The style of local off-fault deformation does not have to be the same as on nearby faults but is usually similar because it is likely to be driven by the same large-scale external forces.We construct a set of synthetic off-fault deformation-style observations (with uncertainties) by regional smoothing of on-fault observations (Figure 3) to guide our fitting process toward such a solution.
We construct a model that fulfills strain compatibility.This is the requirement that strain rates correspond to the derivatives of a velocity field: where the slip rate on a fault dies out, the local velocity field remains compatible through increases in off-fault strain rate.We achieve this by solving for velocities at grid points in the model and calculating strain rates from that velocity field.Two velocities are solved for at each point on a fault, one on each side of the fault, with the fault slip rate being the difference between the two velocities.
Our method combines dynamical principles with those of maximum likelihood regression.Accordingly, we calculate two solution types and iterate until a stable, optimized solution is found.The first step involves computation (by minimization of work rate) of a dynamical forward model (faulted thin viscous sheet) with parameters defining forces and material properties that must be assumed a priori.The second step involves fitting observations to a new model with the same input parameters to make a posteriori inference of an improved kinematic model: a guided non-linear regression that minimizes a sum of squared residuals to observations.By comparing the models, we improve the choice of input parameters until an optimized fit is found.The dynamical model is non-unique, for example, a fault might move faster through changes to either material properties or tractions.We start with an initial guess of material properties with uniform off-fault strength and fault strengths inversely proportional to fault slip rates.We then adjust force potentials to drive faults at the observed rate, and we finally adjust fault strengths to improve the fit.However, the exact choice of formulation is not critical.Our objective in this paper is not to obtain an understanding of the physics, which requires additional observations and/or assumptions.Rather, we construct a dynamical model that follows physical principles, and we then infer a kinematic solution (our primary objective) that departs from physical principles to fit observations, but the kinematic solution must still comply with boundary conditions and strain compatibility.The dynamical model guides the kinematic solution where observations are lacking, that is, interpolation between observations is influenced by the dynamical solution.

Model Boundary Velocities
We specify velocities on the boundaries of our model and require them to fit exactly.The boundaries are set away from the region of interest to minimize potential boundary effects.The model's eastern and western boundaries are set as rigid boundaries where velocities are the rigid body rotations of the Pacific and Australian plates, respectively.We use velocities specified by MORVEL (DeMets et al., 2010).
Velocities on the northern and southern boundaries are specified by interpolating between the two rigid plate-motion rates.Velocities on the northern boundary are set with the Havre Trough having an opening rate of 20 mm/ yr and obliquity defined by an axis of principal extension with azimuth 135° (Caratori Tontini et al., 2019;Wright, 1993).The remainder of the northern boundary is set with no deformation, except for the Hikurangi subduction zone, which is set to accommodate all remaining plate motion.Velocities on the southern boundary are set with Puysegur Ridge having 14 mm/yr dextral strike-slip and the Puysegur subduction zone accommodating all remaining motion.

Faults in the Model
Fault geometries and slip rates used in the model are based on the active fault model of Litchfield et al. (2014).Only larger faults are included in the model as separate faults for which slip rates are modeled (Figure 2).The categorization of larger faults is based on consideration of their slip rate, length, and intersections with other faults.We simplify the geometry of these faults compared to Litchfield et al. and add intersections between faults where they almost intersect (and likely do at depth).Smaller and/or closely spaced faults are combined with neighboring faults and converted to composite strain-rate observations, allowing us to create a manageable and more meaningful (in a thin-sheet sense) model geometry.The accretionary wedge offshore of Fiordland is simplified to a single fault located at the front of the wedge with a shortening rate of 2.8 ± 1.5 mm/yr (Barnes et al., 2002).
We use slip rates given by Litchfield et al. (2014) as our primary set of observations.Litchfield et al. gave rake-parallel slip rates and we convert them to horizontal slip rates using their fault dips and propagated slip rate and dip uncertainties.Each fault segment in the model is assigned a slip-rate observation (observations are halved on fault-terminating segments).For low-sliprate faults where Litchfield et al. (2014) do not assign a slip rate, we assign a nominal slip rate of 0.1 ± 1.0 mm/yr, so that the fault has the correct sense of motion (rake) in the model but its rate is free to vary.Similarly, for faults that do not have a rake assigned, we assign one based on the fault's sense of motion (reverse, strike-slip, etc.) and apply an uncertainty of ±40°.Geological observations show that the Akatore Fault has a time-variable slip rate (anomalous Holocene rate), and we use a slip rate of 0.04 ± 0.01, which is representative of the long-term slip rate on the fault (Taylor-Silva et al., 2020).
The Hikurangi and Puysegur subduction thrusts are large faults with imprecisely determined slip rates.To reflect this, we use the slip-rate estimates of Litchfield et al. (2014) but with larger uncertainties that allow the subduction thrusts to accommodate plate motion that is not easily accommodated elsewhere in the model.For Puysegur, the uncertainty is four times Litchfield et al.'s slip-rate estimate.For northern and central Hikurangi, the uncertainty is half Litchfield et al.'s slip rate estimate and the Hikurangi southern termination (Wellington segment) is left unaltered.
We increase uncertainties by a factor of five from Litchfield et al. (2014) on several high slip-rate offshore faults in the Hikurangi-Marlborough transition: the Needles, Chancet, Campbell Bank, Boo Boo, and Uruti Basin Faults.In addition, we increase uncertainty in the strike-slip component of motion on the Palliser-Kaiwhata Fault by the 7 of 22 same factor, because the reverse-slip component is known from coastal uplift (Berryman et al., 2011), but the strike-slip rate is poorly known.The increase in uncertainty we apply to these offshore faults allows them to vary from their defined best slip-rate value more than others around them, which we suggest is reasonable given that they have abundant geomorphic evidence for a high slip rate, but their true slip rate and uncertainty is unknown.
Large faults and their slip-rate observations are included explicitly in the model, but observations from small and/ or closely spaced faults are included as strain rate observations.We convert slip rates to strain rate using Kostrov summation.The strain rate is where   is the area of the element,   the fault length in that element,  [] fault slip rate and   the normal to the fault (Holt & Haines, 1995).The uncertainties in the slip-rate observations are used to calculate standard errors in the strain-rate observations except where this uncertainty is smaller than a minimum value  2 , in which case the uncertainty is set to  2 .The parameter   , which is varied between 10 −11 /yr and 10 −4 /yr, is explained below.
A uniform strain rate of 1.5 ± 0.4 × 10 −7 /yr is estimated over the width of Havre Trough, to match the 20 ± 5 mm/ yr of oblique extension (Caratori Tontini et al., 2019;Wright, 1993).Regions with negligible strain rates, based on seismic-reflection mapping observations, are assigned observed strain rates of 0 with an uncertainty of 0.02   .These regions are offshore west of New Zealand, offshore east of New Zealand, and the Raukumara Basin (Figure S1 in Supporting Information S1).Other regions without modeled small faults are assigned strain-rate observations of 0 with an uncertainty of  2 .

Style and Direction of Distributed Deformation
In order to guide the style and direction of off-fault deformation, we create a set of strain-rate style observations that we fit simultaneously with on-fault slip observations.Fault slip rate is spread across the region surrounding the fault using a 25 km triangular taper to give an equivalent strain rate.Contributions from all faults (Equation 1) are summed to create a strain-rate field with a value at all points in the model.The reference strain rate in each triangular element is the mean of the strain rates at its vertices.
We then quantify the strain rate style,  0 (Figure 3), which is defined as the dilatation divided by strain-rate magnitude (adapted from Klein et al. [2009]) , representing isotropic expansion, but such circumstances are rare (e.g., around a volcano).In the more common tectonic context, values of  0 from zones of sub-parallel faults range from −1, representing pure reverse faulting, to 1, representing pure normal faulting.
0 = 0 represents pure strike slip faulting.Values outside the range −1 to 1 can occur where there are contributions from faults with different orientations or slip directions.Values of  0 are extrapolated to areas with no faults to provide values of  0 for the entire model domain.The reference values of  0 determined from faults are included as a set of strain-rate observations, but separately from observations of strain rate determined from small and/or closely spaced faults.To implement strain-rate style observations, we use strain-rate magnitude extracted from the previous iteration and multiply it by the reference  0 value to compute an observation of the dilatational component of strain rate  ( ė + ė)  = 0 ėmag .We set the uncertainty in  0 to ±1 for all elements to reflect that regions surrounding a fault may not exactly replicate the style of the fault but are unlikely to follow a different tectonic regime.In our first iteration, the strain-rate magnitude is set as 0.1 times the magnitude of the reference strain rate.Our approach allows strain-rate style to be fit with minimal influence on strain-rate magnitude in the final model solution.
We use the same reference strain-rate field from faults to determine reference directions   of maximum horizontal shortening (Figure 3).To better fit the direction of maximum horizontal shortening, we rotate the strain rate by   , which provides the relationship and then we use Equation 3to construct an observation with uncertainty  Δ  mag .We set the uncertainty in direction  Δ to ±90° for all elements and   mag is the same as used for the strain rate style observation.Equation 3 can be satisfied by both   and  + ∕2 .To distinguish between the two directions, we use as an observation the difference between the maximum and minimum principal strain rates   max and   min : We use an uncertainty of  √ 2 −  2 0 ėmag∕2 to ensure that the observation of Equation 4 has the correct sign.As with strain-rate style, our approach allows strain rate maximum horizontal shortening direction to be fit with minimal influence on strain-rate magnitude in the final solution.

Solution Method
We use codes provided by Haines and Sutherland (2018) to solve for our kinematic model.We present an outline of their method below and provide further details in the appendix, but refer the reader to Haines and Sutherland (2018) for a full description of their method.

Observation Fitting
We seek a kinematic model that fits observations and interpolates realistically between observations.We achieve this through minimization of The ) by dynamical principles and a priori parameters.The approach allows diverse observations with different uncertainties to be included that influence the model in an unbiased way.Each component of an observation (e.g., strike-slip or dip-slip fault slip rate) makes a separate contribution to the weighted least squares statistic, and the approach allows partial observations to be used, for example, just the dilatational component of strain rate or the dip-slip component of fault slip rate.
In some regions, sparse observations may leave the model under-determined and physical principles drive the solution, whereas in other places, it may not be possible to fit all observations and the model is over-determined by observations.This would, for example, occur if the slip-rate observations summed to more than the total plate motion and the observations of deformation style were of the same style as the faults.In order to fit the style observations, the modeled fault slip rates would need to be smaller than their observations.This would create a trade-off between fitting deformation style and fitting slip rate observations.Increasing the weight (i.e., reducing  obs ) of deformation-style observations would increase the model's fit to them but reduce the model's fit to slip-rate observations.

The functional
PRE is the work rate of the physical system and influences how velocities are interpolated between observations: deformation in under-constrained regions is guided by dynamical principles.The functional  PRE HIRSCHBERG AND SUTHERLAND 10.1029/2022JB024828 9 of 22 depends on a priori parameters of material properties and force potentials that are specified at the start of each iteration.It takes the form ( ė) is the work rate related to strain rates    in model elements, with each element receiving a weight   , which is a material property called the strain-rate capacity (inverse of viscosity). ([]) represents work contributions from slip rates  [] on fault segments in the model, with each segment receiving a weight   , which is a material property called the slip-rate capacity (inverse of fault strength; see Appendix).The strain rates and slip rates are affected by the boundary condition, material properties and force potential terms   .
Weighting between statistical and dynamical parts of the solution is achieved through the scale factor   .A larger value of   results in a solution that more closely fits dynamical principles and has a smoother velocity field, at the expense of a poorer fit to observations.A smaller value of   will result in a solution that better fits observations but departs more from dynamical principles and has a velocity field with greater local variations.A reference value   can be defined so that each observation and each degree of freedom in the model are statistically expected to have equal weight (see Appendix).We use an   that is 0.1 times the reference value, which weights by a factor of 10 the fitting of observations over consistency with dynamical principles (physics).It is primarily a kinematic solution that we wish to obtain, but we retain physical plausibility through inclusion of strain-rate style observations and having  = 0.1 .

Dynamical Principles
We use dynamical principles to interpolate velocities between observations, however, as elaborated on below, the final kinematic model is not dynamically consistent.The dynamical relationship between local strain rates    in an element, stress   and strain-rate capacity   is This local relationship is a Newtonian approximation, but every element has its own strain-rate capacity, so the global rheology can be controlled during an adjustment process and may vary spatially.This subject is outside the scope of this paper.We use a uniform strain-rate capacity.models testing values between 10 −11 /yr and 10 −4 /yr (multiplied by units of stress).We also use   to determine the minimum uncertainty in strain-rate observations as larger values of both correspond to more off-fault deformation.

Slip rates
The dynamical a priori parameters of force potentials and material properties are adjusted iteratively using a two-step process.Equations 8-10 are solved through minimization of   to create an a priori dynamical solution.Minimizing  POST causes the model to deviate from the dynamical solution to fit observations better but it maintains strain compatibility and the boundary condition.The difference between the two solutions provides a basis for choosing a new set of values for parameters and hence a new iteration starts.The process is continued until a stable solution is obtained.

Parameter Optimization
The local differences in kinematics between solutions that minimize  POST and  PRE provide an indication for how to adjust each local parameter value, because it reveals the kinematic difference required to better fit observations.Each adjustment of a priori parameters moves the next iteration's a priori solution closer to the preceding iteration's a posteriori solution.Because a posteriori solutions better fit observations than their corresponding a priori solution, this results in a progressively improved fit to observations in both the a priori and a posteriori solutions.
Parameter optimization could be approached through adjustment of capacities   and   or adjustment of forcing potentials   .For example, if a fault is predicted to have a smaller slip rate in the a priori solution than in the a posteriori solution, the slip rate can be increased by increasing the fault slip-rate capacity   or by changing force potentials such that the traction on the fault in the direction of fault slip is increased.As we are seeking a kinematic solution, either of these approaches would be valid, as they both could improve the fit to observations.However, if a fault is predicted to slip in a different direction in the a priori solution than in the a posteriori solution, adjusting   will not suffice as a change in slip direction implies making   negative, which has no physical relevance and hence we do not allow.In this situation, an adjustment to   is needed to drive the fault in the observed direction.This has relevance to, for example, back-arc extension in a subduction zone with overall shortening.Force potentials are required to drive back-arc extension in North Island.
We start by adjusting force potentials   and leave capacities   and   unchanged.The adjustment to   in an element is to add the difference between the values of   ∕ in the a posteriori and a priori solutions.Similarly, the adjustment to   on a fault is the difference between the values of  []∕ in the a posteriori and a priori solutions.The derivation for this relationship is presented in the appendix.This adjustment is mathematically equivalent to adjusting   by the difference between the two solutions in stress   in elements or in traction on faults in dynamical modeling.
To ensure finite forces, force potentials must be smooth.This is achieved in the code by specifying force potentials at points and determining force potentials acting on segments and elements from those points.The force potential adjustment at a point is calculated as the mean adjustment of the surrounding elements.If the point is on a fault, we give it priority and the force potential adjustment is the mean adjustment of surrounding fault segments.
Because force potentials on faults are shared by fault segments and the elements surrounding the faults, adjustments to force potentials based on slip rates also affect strain rates in neighboring elements.The size of the effect on strain rates depends on the relative values of slip-rate capacity and strain-rate capacity, with a larger relative slip-rate capacity resulting in smaller effects on strain rates next to the fault.Adjustments to   allow the size of this effect to vary to further improve the fit to observations.After reaching an optimal solution through adjustments to force potentials, we improve the fit by adjusting slip rate capacities   .The adjustment to   on a fault segment is to multiply it by the ratio of the slip rates in the a posteriori and a priori solutions (see Appendix for derivation).As it is the relative values of   and   that are important rather than their absolute values, it is only necessary to adjust one, so we adjust only   and leave   uniform.The overall solution method is summarized as: 1. Setup: Specification of observations of slip rates, strain rates and strain-rate styles., and observations of slip rates, strain rates and strain-rate styles to create an a posteriori inference of a velocity field that fits observations, guided by but not exactly satisfying physical principles.4. Adjustment of   : Adjustment of force potentials   based on differences between the a posteriori solution and a priori solution.5. Calibration of strain-rate style observations: Magnitude of strain-rate style observations recalibrated to use strain-rate magnitudes from the previous a posteriori solution.6. Iteration: Iteration of steps 2-5 until the a posteriori solution stabilizes at an optimal solution.7.

Comparison of Models
We construct models with   varying from 10 −11 /yr to 10 −4 /yr and optimize to find the a posteriori quality of fit to observations of fault slip rate   (see Figure 4).  is large at ∼30/observation for   = 10 −11 /yr but decreases rapidly to ∼5/observation for   between 10 −10 /yr and 10 −8 /yr.  decreases smoothly to reach a minimum of ∼0.3/observation at   = 2 × 10 −6 /yr, and then increases again to reach ∼50/observation at   = 10 −4 /yr.Models with small   have relatively stiff elements between faults and, therefore, tend to predict a smaller proportion of plate motion being accommodated off faults (i.e., smaller rates of distributed strain).As the models are constrained by boundary conditions, these models have a tendency to over-estimate fault slip rates.Stiffer elements are also less able to accommodate rapid spatial variations in strain rate between faults, meaning that some observations cannot be fit while satisfying strain compatibility, resulting in a poorer fit to observations and larger   .
As   increases, the elements become weaker and better able to accommodate "missing" plate displacement not accommodated on faults.Weaker elements are also better able to accommodate the strain rate around faults where the slip rate changes over a short distance.This results in a general decrease in   and improvement in fit to fault observations.Increasing   also increases the local contrast between slip-rate capacity and strain-rate capacity.Force potentials at points on faults affect the slip rate on the fault and the strain rate in elements next to the fault in proportion to their respective capacities.At large   , the force potential required to cause a small change in fault slip rate results in a large change in element strain rate.At very large   , this effect produces very large strain rates that are not dampened by other constraints in the model and results in a larger   and poorer fit to observations.This becomes noticeable for   ≥ 5 × 10 −6 /yr.The best-fitting model occurs for   = 2 × 10 −6 /yr and we describe the results of this model, which we call HS22-Lit.In our discussion, we compare it to the model with   = 10 −7 /yr (see Figure S3 in Supporting Information S1).This second model has   of ∼1.1/observation, indicating that the observations are fit at about their expected level of uncertainty.

Overview of Best-Fitting Model
The method described above produces a kinematic model that fits observations (Figures 5 and 6) and the major features of the New Zealand plate boundary (Figure 7).The model balances fault slip rates within their uncertainties and ensures that modeled slip rates are consistent with surrounding fault slip rates (Table 1) and strain rates (Figure 8).Furthermore, by balancing the plate-motion budget, the model uses the combined observations to improve individual inferences of slip rate (and uncertainty) on faults (Table S1), including the subduction thrusts.Subduction thrust slip-rate observations currently have large uncertainties, so improving inferences of slip rates improves inferences of seismic hazard.The model confirms, as expected, that the majority of plate motion is accommodated on the Hikurangi subduction zone, Marlborough fault system, Alpine Fault, and Puysegur subduction zone.There is extension in western and central North Island, shear along NIDFB, MFS, and Alpine Fault, and shortening for most other parts of the boundary (Figure 9).Strain rates are generally <0.05 × 10 −6 /yr, but there are larger strain rates in central Island and offshore eastern North Island.

Hikurangi Subduction Zone
The HS22-Lit model has a southwards decrease in the slip rate of the Hikurangi Subduction thrust.Our revised prediction of shortening rate decreases from 33 ± 6 mm/yr in the northern Hikurangi margin to 13 ± 6 mm/yr in the center.
Including deformation near the base of the accretionary wedge, we model a total of 45 ± 8 mm/yr of shortening in northern Hikurangi and 34 ± 3 mm/yr in central Hikurangi.Our model also predicts significant dextral deformation of 15 ± 3 mm/yr in the north and 19 ± 5 mm/yr in the center occurring in the accretionary wedge (though see discussion below about comparison to GPS velocities).In the southern Hikurangi margin, we do not alter the uncertainty from Litchfield et al. (2014) because the smaller uncertainty is required to ensure a reverse component of slip in the model.For this segment, we model 22 ± 5 mm/yr of shortening.Including deformation in the accretionary wedge, we model 20 ± 3 mm/yr of shortening and 23 ± 5 mm/yr of dextral slip in southern Hikurangi.

Wellington-Marlborough
In the transition from the Hikurangi subduction zone to the MFS, slip from the subduction zone is transferred onto a series of offshore faults.The east-west striking Boo Boo Fault in the HS22-Lit model has 10 ± 2 mm/yr of dextral slip.The Needles Fault is modeled with 14 ± 4 mm/yr of dextral slip and 2 ± 3 mm/yr of shortening.We model 18 ± 4 mm/yr of dextral slip and 1 ± 3 mm/yr of shortening on the offshore segment of Kekerengu Fault, from where slip is transferred from onshore faults.The relatively small model standard errors for offshore fault slip rates (with very poorly known slip rates and uncertainties) emerge as a result of strain compatibility with onshore faults and the boundary condition (plate-boundary budget).
In Marlborough, the Wairau Fault in the HS22-Lit model has 3.9 ± 1.1 mm/yr of dextral slip.We model dextral slip of 5.3 ± 1.7 mm/yr on Awatere Southwest and 5 ± 2 mm/yr on Awatere Northeast.The Clarence Fault is modeled with 3.0 ± 0.8 mm/yr of dextral slip for the southwest segment and 4.0 ± 0.7 mm/yr for the northeast segment.The onshore segment of the Kekerengu Fault is modeled with 14 ± 4 mm/yr of dextral slip and 6 ± 4 mm/yr of shortening.Jordan Thrust is modeled with 7 ± 5 mm/yr dextral and 14 ± 3 mm/yr of shortening.The Kelly Fault is modeled with 18 ± 4 mm/yr of dextral slip, transferred as 6.5 ± 1.9 mm/yr on the Kakapo Fault and 13 ± 2 mm/yr on the western portion of the 1888 segment of the Hope Fault.Farther east, the Conway section of the Hope Fault is modeled with 22 ± 3 mm/yr of dextral motion and 2 ± 2 mm/yr of shortening.

Alpine Fault to Puysegur
On the northern sections of the Alpine Fault, where most plate motion is accommodated in the MFS, the dextral slip rate increases from 5 ± 3 mm/yr where it meets Wairau Fault in the north to 18 ± 5 mm/yr where it meets the Kelly Fault in the south.The Alpine Fault accommodates ∼80% of the plate motion in its central and southern sections.We model 29 ± 5 mm/yr of dextral slip and 2 ± 4 mm/yr of shortening for the central section.The dextral slip rate decreases to 23 ± 2 mm/yr just north of Fiordland, but increases again to 31 ± 3 mm/yr adjacent to southern Fiordland.
The Puysegur subduction margin is modeled to have a dextral slip rate of 22 ± 15 mm/yr and a shortening rate of 26 ± 11 mm/yr.The Puysegur Ridge Fault is modeled with a dextral slip rate of 10 ± 11 mm/yr and the Snares Fault is modeled with a dextral slip rate of 11 ± 6 mm/yr.The Fiordland central wedge is modeled with 4.2 ± 1.3 mm/yr of shortening.

Uncertainty in Results
We use Monte Carlo simulations of Litchfield et al.'s (2014) uncertainties to estimate uncertainties in our HS22-Lit model.Our uncertainties are affected by the geometry of the model, which requires strain compatibility (fault slip rates and local strain rate must add up) and consistency with plate motions.Central Limit behavior during the regression results in a more Gaussian distribution of our output uncertainties, which gives our uncertainties greater statistical meaning than the uncertainties of Litchfield et al.Our simulations also allow an estimate of correlations between fault slip rates that can be used in conjunction with our slip-rate estimates.
We estimate uncertainties in modeled velocities (relative to the Pacific or Australian plate) to be ∼5 mm/yr for most of the more active part of the plate boundary, covering the Havre Trough, NIDFB, MFS, and Alpine Fault.Larger uncertainties of ∼10 mm/yr are estimated for the Puysegur and Hikurangi subduction zones, where there We select the HS22-Lit model with our preferred value of   as the model with the best fit to fault slip-rate observations.Velocity differences between this model and the model with   = 10 −7 /yr (Figure S3 in Supporting Information S1), which does not overfit observations, are generally small compared to the uncertainty in the models.This indicates that the precise choice of initial input parameters does not greatly affect the modeled velocities at the scale of tens of kilometers, however, it does affect whether the deformation is modeled as slip rate on faults or strain rate off faults.

Comparison to GPS
The long-term velocity field, as determined in our best-fitting model, displays some notable differences from the contemporary velocity field determined from GPS observations (Beavan et al., 2016;Haines & Wallace, 2020; see Figure 10).GPS observations are not available offshore, so we do not interpret differences offshore where contemporary velocities are poorly constrained.
In southwestern and central western South Island, east of the Alpine Fault, long-term velocities are more southwestward by up to ∼15 mm/yr, with larger differences occurring closer to the Alpine Fault.Long-term rates are closer to the Pacific plate rate and indicate a more localized zone of deformation near the Alpine Fault than in the contemporary GPS.This may be interpreted as a locking signal on the southeast-dipping Alpine Fault (Haines & Wallace, 2020;Lamb & Smith, 2013;Wallace et al., 2007).While some of the difference may be due to deformation not occurring on faults in our HS22-Lit model, we only predict ∼9 mm/yr of long-term deformation (see Section 5.4) so this cannot explain all differences.In southern North Island, contemporary velocities have a west component larger by up to ∼20 mm/yr, making them closer to Pacific plate velocity than long-term HS22-Lit velocities.This indicates that contemporary movement in this region includes a component derived from the Hikurangi subduction zone and that the southern portion of the Hikurangi subduction zone is strongly coupled, consistent with previous suggestions (e.g., Walcott, 1984;Wallace et al., 2004Wallace et al., , 2012)).Strong coupling has been interpreted as corresponding to a slip deficit, most of which may eventually be released through an earthquake or aseismic slip on the subduction interface (e.g., Jolivet et al., 2020;Moreno et al., 2010;Wallace et al., 2004).
In eastern North Island and the northern Hikurangi subduction zone, contemporary velocities have a larger south component by up to ∼10 mm/yr.Our HS22-Lit model long-term rotation of eastern North Island is smaller than observed in contemporary velocities.Rotation of northeastern North Island is required by the southward decrease of extension rates in the TVZ and the Havre Trough (Beanland & Haines, 1998;Wallace et al., 2004).Two possibilities arise for why our predicted rotation of northeastern North Island is less than geodetic estimates: (a) the geodetic signal is transient and future large earthquakes will result in a long-term velocity field that is different to the short-term velocity field; or (b) our long-term velocity field is inaccurate.We favor the second interpretation, because it seems unlikely that the transient dextral strike-slip signal (for north-striking faults) is consistent with subduction zone locking.Instead, we suggest that our velocity field may be inaccurate because dextral slip in Bay of Plenty and NIDFB may be partly obscured by poor exposure and/or because there is high uncertainty in the rate and direction of extension in Havre Trough.These suggestions require testing by detailed field investigations, because 5-10 mm/yr of missing dextral fault motion in the central North Island has significant implications for local seismic hazard in a region with moderate vulnerability.

Hikurangi Subduction Zone
Slip rates for the Hikurangi and Puysegur subduction thrusts have previously been estimated from regional kinematic constraints, but have high uncertainty attached (e.g., Litchfield et al., 2014;Wallace et al., 2004).The large uncertainties we input for subduction margin slip rates allow our model to determine best-fitting slip rates based on the plate motion that is not accommodated in other ways.Including shortening in the offshore accretionary wedge, our HS22-Lit model has 45 ± 8 mm/yr of shortening in the northern margin and 34 ± 3 mm/yr in the center.For both parts of the margin, this is smaller than the respective values of 54 ± 6 mm/yr and 44 ± 7 mm/yr modeled by Wallace et al. (2004), but we identify the same pattern of shortening decreasing southwards.Additionally, we model significant dextral motion of 15 ± 3 mm/yr in the north and 10 ± 6 mm/yr in the center.We model this as being mostly accommodated in the lower trench slope, but it may be accommodated further up the trench slope or within the NIDFB.

Our model with
= 10 −7 /yr predicts larger shortening rates of 36.7 ± 1.3 mm/yr and 23.6 ± 1.2 mm/yr on the northern and central portions of the Hikurangi subduction thrust.Note that the standard errors on the prediction are smaller, because the effective number of degrees of freedom is higher in a model with smaller   , because distributed deformation is more constrained by higher material strength.The model with   = 10 −7 /yr predicts larger dextral slip rates of 7.9 ± 1.3 mm/yr and 11 ± 2 mm/yr on the northern and central subduction thrust, respectively.Including contraction in the accretionary wedge, this model predicts 42.1 ± 0.7 mm/yr of shortening and 15 ± 2 mm/yr of dextral slip in the northern portion; and 31.6 ± 0.9 mm/yr of shortening and 18.2 ± 1.4 mm/yr of dextral slip in the central portion.The rates for this part of the subduction zone, which includes the accretionary wedge, are not significantly different to those predicted by the best-fitting model but have smaller uncertainties due to the higher number of effective degrees of freedom.While our HS22-Lit model has 22 ± 5 mm/yr of shortening on the southern portion of the Hikurangi subduction thrust, we model a smaller amount of 20 ± 3 mm/yr of shortening across the entire accretionary wedge.To accommodate this difference, some extension is modeled in the accretionary wedge, indicating that observations nearby add up to more shortening than is available in the plate boundary budget.The model with   = 10 −7 /yr models a notably smaller shortening rate of 4 ± 2 mm/yr on the subduction thrust, along with 11 ± 4 mm/yr of dextral slip.Including deformation in the accretionary wedge, this model predicts a smaller shortening rate of 10.5 ± 0.9 mm/yr but its dextral deformation rate of 23.8 ± 1.2 mm/yr is not significantly different from the best-fitting model.

Alpine Fault
The intersection between the Alpine, Kelly and Hope Faults is a region of relatively large misfits in the HS22-Lit model.The Kelly Fault is modeled with a dextral slip rate of 18 ± 4 mm/yr compared to 13 ± 2 mm/yr estimated by Litchfield et al. (2014).Between the Kelly and Hope Fault intersections, the Alpine Fault is modeled with 18 ± 5 mm/yr of dextral slip, decreasing to 15 ± 3 mm/yr north of the Hope Fault intersection, larger than the 11 ± 2 mm/yr specified by Litchfield et al.This indicates that the estimated fault slip rates in this region are inconsistent with each other, so estimated rates are unable to match all observations (because of strain compatibility).Further observational work there is required.
The HS22-Lit model dextral slip rates on the central Alpine Fault of 29 ± 5 mm/yr are consistent with the observed slip rate of 28 ± 4 mm/yr (Howarth et al., 2018;Norris & Cooper, 2001).Including strain rate modeled adjacent to the Alpine Fault, 38 ± 7 mm/yr of motion is modeled parallel to the Alpine Fault.The model with   = 10 −7 /yr predicts a larger slip rate of 34.4 ± 1.4 mm/yr on central Alpine Fault but predicts small strain rates adjacent to the fault.Overall, it predicts a total of 38.9 ± 0.3 mm/yr of dextral strike slip on or near (within ∼50 km) the Alpine Fault, consistent with the prediction of the best-fitting model.Southeast of the Alpine Fault are several faults in the Southern Alps with unknown slip rates due to rapid erosion occurring in the mountainous region (Cox et al., 2012).Furthermore, the fault model may be incomplete in this region (Litchfield et al., 2014), meaning the ∼9 mm/yr of motion as distributed strain rate in the HS22-Lit model may be accommodated by a combination of known faults with unknown slip rates, currently unknown faults, and/ or off-fault deformation.

Puysegur-Fiordland
The Puysegur subduction thrust has a HS22-Lit model dextral slip rate of 22 ± 15 mm/yr and a shortening rate of 26 ± 11 mm/yr.While these rates are consistent with the strike slip rate of 11 ± 4 mm/yr and shortening rate of 23 ± 8 mm/yr specified by Litchfield et al. (2014), our model predicts that the subduction thrust has more oblique slip and accommodates a greater portion of the plate motion.The Puysegur Ridge Fault is modeled with a dextral slip rate of 10 ± 11 mm/yr compared to the 14 ± 14 mm/yr specified by Litchfield et al.The Snares Fault is modeled with a dextral slip rate of 11 ± 6 mm/yr compared to 10 ± 10 specified by Litchfield et al.Slip rates in this region are poorly constrained and the uncertainties are typically 0.5-1 times the slip rate magnitude.

The New Zealand Community Fault Model v1.0
We repeat our modeling using the recently released New Zealand Community Fault Model (NZ CFM) v1.0 (Seebeck et al., 2022;Van Dissen et al., 2021) in place of Litchfield et al.'s (2014) model.The NZ CFM-based model, which we call HS22-CFM, has a mean fit to fault observations of ∼0.8/observation.The modeled slip rates are provided in Table S2.HS22-CFM is similar to our original HS22-Lit: velocity differences are typically <5 mm/yr (Figure S4 in Supporting Information S1).This similarity is expected as the NZ CFM draws heavily on the Litchfield model.The largest differences between HS22-CFM and HS22-Lit occur in eastern North Island and the Hikurangi trench slope where HS22-CFM predicts larger eastward velocities.HS22-CFM shortening rates in the wedge are modeled at 44 ± 4 mm/yr in the north and 33 ± 4 mm/yr in the center, which are consistent with rates modeled in HS22-Lit.The shortening rate of 10 ± 2 mm/yr modeled in the south is approximately half the rate in HS22-Lit but is consistent with the rate from the model with   = 10 −7 /yr.The dextral deformation rate of 8 ± 3 mm/yr modeled in the north is approximately half the rate in HS22-Lit whereas the dextral rate of 17 ± 3 mm/yr in the center is larger than, but still consistent with, HS22-Lit.The dextral rate of 25 ± 3 mm/yr in the south is consistent with HS22-Lit.
HS22-CFM predicts velocities slightly closer to contemporary velocities in eastern North Island than HS22-Lit does, but still with differences of up to ∼8 mm/yr for HS22-CFM.HS22-CFM still predicts velocities that are less southward than contemporary velocities.The uncertainty in modeled velocities in eastern North Island is ∼5 mm/yr.
In HS22-CFM, we do not alter the Akatore Fault from the faster, Holocene horizontal slip rate of 0.9 (0.2-4.5) mm/yr given in NZ CFM v1.0.However, HS22-CFM models a significantly smaller horizontal slip rate of 0.05 mm/yr, consistent with the smaller, longer-term slip rate of 0.04 ± 0.01 (Taylor-Silva et al., 2020) that we used in HS22-Lit.This implies that kinematic and strain compatibility constraints likely require the Akatore Fault to average the slower rate over longer timescales.

Conclusions
We have presented a model of deformation across the New Zealand plate-boundary zone that fits Quaternary fault slip-rate observations and is constrained by known plate boundary motion.Our model includes off-fault deformation that fits with our model of deformation style derived from nearby faults.Unlike previous models, we include discontinuities at faults in the model, allowing us to improve the resolution of the model from ∼100 km to <10 km and hence to estimate subduction zone slip rates.
Strain compatibility allows us to identify regions where current models of slip rate are incompatible.One such region is the intersection between the Alpine, Kelly and Hope Faults, where comparatively large misfits suggest that the current models of slip rate should be re-evaluated there.Our HS22-Lit model indicates that there is ∼9 mm/yr of plate motion in central South Island not accounted for by current models of fault slip rates.This may partially be accommodated by faults with unknown slip rates immediately southeast of the Alpine Fault.We model deformation in this region as occurring over a narrower band than indicated by contemporary GPS-derived velocities, which is poorly constrained by observations but expected for elastic loading on a southeast-dipping fault.
Our HS22-Lit model predicts shortening rates across the Hikurangi subduction margin that are significantly smaller than previously published estimates.In the southern Hikurangi subduction zone differences to contemporary velocities are consistent with a strongly coupled subduction thrust and it has previously been suggested that this deficit will be reconciled by a future slip event on the southern Hikurangi subduction thrust.In contrast, the contemporary velocity field in northeast North Island is 5-10 mm/yr southward relative to our modeled long-term velocity field.Based on the mismatch with what is expected for a subduction locking signal, we suggest that this difference arises through inaccuracy of our long-term velocity field, which is potentially due to uncertainty in extension in Havre Trough and/or an underestimate of strike-slip motion in northeastern North Island.
Our analysis provides a quantitative basis for evaluating how much deformation is missing from existing active fault models, and where it might be required.Our results provide a focus for re-evaluating existing field data and will help target collection of new data.Perhaps most importantly, the subduction thrusts in New Zealand are the most rapidly slipping and hazardous faults, but notoriously difficult to quantify; and our method allows their slip rates to be determined from a rigorous analysis of the plate boundary slip-rate deficit.Finally, our method could be applied to any other region where plate-motion boundary constraints are known and there are rich quantitative observations that can be directly related to any component of fault slip rate or crustal strain rate.
Appendix A: Details on Solution Method

A1. Observation Fitting
The contribution based on dynamical principles to the solution is given in Equation 7 in terms of  ( ė) and   ([𝑢𝑢]) .In full, these terms are

A2. Parameter Optimization
The parameters   ,   , and   can be adjusted to better fit observations, with the difference between the a priori and a posteriori solutions providing an indication of an efficient optimization strategy.To demonstrate our strategy, we consider the local effects of small changes in these parameters on the modeled strain rate and slip rate.

Figure 1 .
Figure 1.Map of New Zealand's tectonic setting.Arrows indicate the velocity of the Pacific (PAC) plate relative to the Australian (AUS) plate.Faults used in the model directly are indicated by thick red lines.Faults included in the model as strain rate are indicated by thin purple lines.TVZ, Taupo Volcanic Zone; NIDFB, North Island Dextral Fault Belt; MFS, Marlborough Fault System; Fiord., Fiordland; SZ, subduction zone.

Figure 2 .
Figure 2. Slip-rate and strain-rate observations used in the model.Faults are colored by their slip rate from Litchfield et al. (2014).Elements are colored by the magnitude of their strain-rate observation.Small faults are included in the model by converting their slip rates into strain-rate observations according to Equation 1.

Figure 3 .
Figure 3. Model of strain-rate style and direction used in the model.Slip rates from all faults are converted to strain rates and spread according to Equation 1. Strain-rate style (  0 ) is defined by Equation 2 and extrapolated to points away from faults so that the model of strain-rate style covers the entire model domain.Values of  0 can range from 1.4 representing isotropic expansion (red) to −1.4 representing isotropic contraction (dark blue).Yellow indicates normal (Nor.)faulting and is present in western North Island.White indicates strike-slip (SS) faulting and is present along the Alpine Fault, Marlborough Fault System, and North Island Dextral Fault Belt.Blue indicates reverse (Rev.)faulting and is present along the Puysegur and Hikurangi subduction zones and in Otago.Dark bars indicate the maximum horizontal shortening direction.Light gray indicates regions explicitly made more rigid in the model.

Figure 5 .
Figure 5. Misfit of model HS22-Lit to slip-rate and strain-rate observations.Misfits to slip-rate observations are plotted as blue circles with the diameter of the circle proportional to the misfit.The misfit is plotted in the central segment of the faults, as indicated by the small white circles.Gray squares indicate subduction zones where misfit is considered separately (see main text).Misfits to strain-rate observations are indicated by the colored background.

Figure 6 .
Figure 6.(a) Misfit of model HS22-Lit to strain-rate style, defined as the absolute difference between observed A 0 shown in Figure 3 and modeled A 0 from Figure 9.The model fits the style very well in almost all regions, except for the southern Hikurangi margin.(b) Misfit to strain-rate direction, defined as the absolute difference between the observed direction of strain rate shown in Figure3and the modeled direction from Figure9.The model fits the direction well in most places, but with several localized areas that fit poorly.

Figure 7 .
Figure 7. Modeled velocity fields of model HS22-Lit relative to Australian (a) and Pacific (b) plates demonstrating prominent features of the plate boundary zone.This includes the clockwise rotation of the eastern North Island driven by extension in the Taupo Volcanic Zone and the Havre Trough.The model also shows the Hikurangi and Puysegur subduction zones and Alpine Fault accommodating most of the relative plate motion.Additionally, the effect of dextral faults in the Marlborough Fault System and North Island Dextral Fault Belt is shown.The orange lines indicate locations where we calculate deformation across the Hikurangi subduction margin, including the accretionary wedge.

Figure 8 .
Figure 8. Modeled HS22-Lit strain-rate magnitude.Large strain rates are modeled for extension in the Taupo Volcanic Zone and Havre Trough.Large rates are also modeled for compression in the Hikurangi subduction margin and for strike-slip adjacent to the Alpine Fault.Strain rates of <0.05 × 10 −6 / yr are generally modeled elsewhere.The orange lines indicate the locations where we calculate modeled deformation across the Hikurangi subduction margin, including the accretionary wedge.

Figure 9 .
Figure 9. Modeled HS22-Lit strain-rate style and direction.Style is defined by Equation 2. Yellow indicates extension and blue compression.Bars indicate direction of maximum horizontal compression.Extension is modeled for the Havre Trough and western North Island.Compression is modeled for the Hikurangi and Puysegur subduction zones and Otago.A band of strike-slip is modeled along the Alpine Fault, MFS and NIDFB.Light gray indicates regions that are relatively rigid in the model.

Figure 10 .
Figure 10.Comparison between long-term velocity field from the HS22-Lit model and contemporary velocity field from GPS (Haines & Wallace, 2020).The difference contemporary minus long-term is indicated by arrows and colors show the magnitude of the difference.
Litchfield et al. (2014)s and observations (Ob.) are fromLitchfield et al. (2014), with uncertainties modified on the subduction zones and the other faults noted in the main text.Long.and Lat. are the longitude and latitude of the center of the segment.

Table 1
Summary of Slip-Rate Observations and HS22-Lit Predictions on Selected Fault Segments accommodated off faults, with larger strain-rate variations and smaller slip-rate variations.Conversely, smaller values of   relative to   result in more deformation being accommodated on faults, with larger slip-rate variations and smaller strain-rate variations.Consequently, there is a trade-off where a larger value of   results in smoother slip rates that are a better fit to slip-rate observations (our primary observations) at the expense of a less smooth strain-rate distribution.We use a simple formulation with   initially set in proportion to the slip-rate observation on fault segments (multiplied by units of stress).We set a uniform Iteration with adjustment of   : Repeat steps 2-6, replacing adjustment of   in step 4 with adjustment of slip-rate capacities   based on differences between a posteriori and a priori solutions.
( ė) = ∫ − ∑ ( ė)represents the contribution from strain rates in elements and is integrated over the surface of the model, excluding faults. ([]) represents the contribution from slip rates on fault segments and is integrated over the lengths of all fault segments.  and   represent weights applied to strain rate    in elements and slip rate  [] on fault segments, respectively.   represents a forcing term that is used to better fit observations.  is the horizontal normal to the fault.The dynamics-based contribution is weighted by the factor   in Equation 5 and is specified in terms of a reference value.The reference value faults∫      ([] − ) ([] − )  (A1b) −  represents the effects of the boundary condition.Neglecting the second order effects in Note that the slip rate only specifies adjustments to the two components of   corresponding to the tractions on the fault.The third component of   does not affect the modeled slip rate and therefore adjustments to it are not constrained by desired changes in slip rate.