Robust Imaging of Fault Slip Rates in the Walker Lane and Western Great Basin From GPS Data Using a Multi‐Block Model Approach

The Walker Lane (WL) in the western Great Basin (GB) is an active plate boundary system accommodating 10%–20% of the relative tectonic motion between the Pacific and North American plates. Its neotectonic framework is structurally complex, having hundreds of faults with various strikes, rakes, and crustal blocks with vertical axis rotation. Faults slip rates are key parameters needed to quantify seismic hazard in such tectonically active plate boundaries but modeling them in complex areas like the WL and GB is challenging. We present a new modeling strategy for estimating fault slip rates in complex zones of active crustal deformation using data from GPS networks. The technique does not rely on prior estimates of slip rates from geologic studies, and only uses data on the surface trace location, dip, and rake. The iterative framework generates large numbers of block models algorithmically from the fault database to obtain many estimates of slip rates for each fault. This reduces bias from subjective choices about how discontinuous faults connect and interact to accommodate strain. Each model iteration differs slightly in block boundary configuration, but all models honor geodetic and fault data, regularization, and are kinematically self‐consistent. The approach provides several advantages over bespoke models, including insensitivity to outlier data, realistic uncertainties, explicit mapping of off‐fault deformation, and slip rates that are more objective and independent of geologic slip rates. Comparisons to the U.S. National Seismic Hazard Model indicate that ∼80% of our geodetic slip rates agree with their geologic slip rates to within uncertainties.


Introduction
In the Walker Lane (WL) and western Great Basin (GB) east of the Sierra Nevada Mountains active faults accommodate ∼20% of the Pacific/North America plate boundary relative motion (Bennett et al., 2003;Dokka & Travis, 1990;Oldow, 2003;Thatcher et al., 1999).The fault system is composed of a complex set of active dextral, sinistral, and normal faults (Faulds & Henry, 2008;Stewart, 1988;Wesnousky, 2005) that together work to release the accumulating crustal strain (Figure 1).The faults are numerous, in places closely spaced and discontinuous, and are linked to regional seismicity (dePolo, 2008;dePolo & dePolo, 2012).The system's complexity creates a challenge for using measurements of active deformation from GPS geodesy to estimate the slip rates on the faults.
Fault slip rates are important for several reasons.For example, rates and azimuths of fault slip are used to understand the kinematics of tectonics that drive processes within the plate boundary zones (e.g., Weldon & Humphreys, 1986), and are used to compare to results of global plate tectonic circuits (e.g., deMets and Merkouriev, 2016).Another reason that is of high societal relevance is that slip rates are used in seismic hazard analysis, such as the US National Seismic Hazard Maps (NSHM) project.In these hazard models the faults represent sources of elastic moment that will be released in potentially damaging future ground shaking from earthquakes (Field et al., 2023;Frankel et al., 2002;Peterson et al., 2014).For all these applications the accuracy of the faults slip rates is essential.
However, the problem is challenging owing to gaps in the fundamental data sets, and aleatory and epistemic uncertainties in those data.Geodetic data can constrain the rate of interseismic strain accumulation expected to be released in earthquakes but is limited by the short time of observation (usually one to three decades).Nonetheless, they are a vital complement to other data such as seismicity and earthquake geology which also constrain hazard (Bird, 2009;Kreemer & Young, 2022;Shen et al., 2007).These data types each have their own strengths, sources of uncertainty, community of practice, and offer complementary means for measuring active crustal deformation.Comparisons between geologic, geodetic, and seismic moment rates reveal both similarities between, and gaps among, the data sets (Pancha et al., 2006;Ward, 1998), suggesting that approaches that take advantage of their individual strengths are needed.
Integrated approaches have been attempted in a recent modeling exercise to support the latest version of the US NSHM.A group of expert modelers developed solutions from the same input database of western US fault geometries and GPS velocity field.They estimated slip rates for faults in the database using their own various methodologies incorporating combinations of geologic and geodetic data (Evans, 2022;Pollitz, 2022;Shen & Bird, 2022;Zeng, 2022a).Introduction and review of these models have been provided by Pollitz et al. (2022b) and Johnson et al. (2024).They point out that while the slip rate estimates exhibit broad similarity, they differ in some important measures owing to multiple factors that sometimes result in high variability.In those models the coefficient of variation for slip rates below 5 mm/yr (a category inside which all WL faults lie) is above 2.0, indicating lack of agreement among the modelers at the level of the uncertainties.Moreover, all the models incorporate geologic slip rates (Hatem et al., 2022b) as additional constraints to regularize their inversions and so the resulting estimates are not independent of geologic rates.
Here we address the slip rate problem by introducing a new block modeling method where block geometries are repeatedly generated through a well-defined computational procedure from a given input fault data base.In each  iteration the model starts with many blocks whose number are iteratively reduced to limit the number of free parameters.These models are generated automatically in subdomains of the region of interest, iterating with slightly different starting conditions so that large numbers of slip rate estimates are made for each fault.While each individual model contains some of the usual errors associated with model construction, the errors are different for each iteration and so are a noise that is reduced by averaging over many models.However, the signals of fault slip rates are similar in every model because they are present and have the same geometry in each iteration.In this way the power of large numbers increases the robustness of the slip rate estimate to find the set of kinematically consistent slip rate estimates that fit the data.The technique is like other robust approaches that involve repeated sampling of the data to achieve estimates of velocity or strain rates (e.g., Blewitt et al., 2016;Husson et al., 2018;Kreemer et al., 2018).
The degree of independence that the estimated geodetic rates have from geologic rates is an issue with the NSHM deformation models identified in recommendations for future efforts (Johnson et al., 2024).In our method described here geologic data from the NSHM (Hatem et al., 2022a) are used to define the location, geometry, and style of fault segments, but not their slip rates.Thus, agreements between our geodetic and geologic rates are corroborative, rather than the result of an analytical constraint that they must be similar.

Velocity Data
We use horizontal position time series from the GPS holdings of the Nevada Geodetic Laboratory (NGL, Blewitt et al., 2018) which are processed uniformly using the GipsyX software (Bertiger et al., 2020).RINEX data are from networks listed in the Data Availability Statement and include stations whose time series have duration longer than 2.5 years.We used data products including daily reference frame alignment, orbit and clock files that were provided by the Jet Propulsion Laboratory.More properties of the GPS data processing may be found in Kreemer et al. (2020) and in the NGL GPS data analysis strategy and products summary (http://geodesy.unr.edu/gps/ngl.acn.txt).We correct the time series for the effects of non-tidal atmospheric, non-tidal ocean, and hydrological surface mass loading.These corrections have been shown to reduce noise in GPS time series, especially in the vertical component, but also in the horizontal component (Chanard et al., 2018;Martens et al., 2020).The corrections are based on the predictions from the Earth System Modeling Group of GFZ Potsdam (Dill & Dobslaw, 2013), which are provided on a 0.5°× 0.5°grid and interpolated onto the GPS station locations and provided on the NGL GPS station pages.The trends of the time series represent motion with respect to a fixed North American (NA) plate that has an Euler pole of rotation in the International Terrestrial Reference Frame (Altamimi et al., 2016) defined by Kreemer et al. (2014).We obtain trends in station positions from the corrected time series with the MIDAS robust non-parametric estimator, which is insensitive to steps, seasonality, outliers and heteroskedasticity in the time series (Blewitt et al., 2016).
Additional campaign GPS data are taken from several sources that are listed in the Data Acknowledgement section.While campaign GPS velocities tend to be less precise than those from continuous stations, they have been shown in many studies to be precise enough to constrain crustal movement.They are numerous and dense, enhancing geographic coverage of the velocity field (e.g., Bennett et al., 1996;Hammond & Thatcher, 2005;Lifton et al., 2013;McCaffrey et al., 2007;Murray & Svarc, 2017;Oldow, 2003;Spinler et al., 2010;Thatcher et al., 1999, to name a few).We align these velocity solutions into the same North America reference frame as the NGL MIDAS NA velocity field.Stations in clusters within 0.001°of one another were combined into single rates by taking the median rate of the cluster.We include all stations within 1°outside the model domain (Figure 1) to help constrain the reference frame alignments between velocity fields and to reduce edge effects in the analysis that follows.Together these networks provide 1,311 stations inside the domain (421 continuous, 343 MAGNET semicontinuous, 547 campaign) (Figure 1).Like other recent analyses we impose a velocity uncertainty floor of 0.1 mm/yr for both horizontal components to prevent very low uncertainty velocities from having a disproportionately large influence on the inversions for slip rates.Histograms of the east and north component velocities and their uncertainties are provided in Figures S2 in Supporting Information S1 and in map view in Figure 1b and Figure S3 in Supporting Information S1.A table of velocities and uncertainties is provided in the Supplement (Table S1).
To omit outliers and noise from the velocity field we apply median spatial filtering to each horizontal component.This process replaces each velocity with a weighted median of itself and its neighbors, where the weight is a function of distance, data uncertainties, and a spatial structure function (Hammond et al., 2016).The same process is used to image the velocities by estimating at every point on a regular grid (spacing of 0.045°× 0.045°or ∼5 × 5 km) the weighted median of velocities from the nearest stations.The gridded field is used for the block modeling because it is representative of the smoothly varying WL velocities and ensures that multiple velocity estimates are present for even small blocks.The gridded field is shown in Figure S3 in Supporting Information S1.
Six profiles of velocity across the WL normal to the direction of Pacific and North America relative plate motion show the patterns of the velocities and the magnitude of the signals with respect to their uncertainties (Figure 2).The location of the 50 km wide profiles is shown in Figure 1a.The gridded velocity is always within the uncertainty bounds of the station observations, except in cases with obvious outliers.The profile-perpendicular velocities increase to the west monotonically and their spatial gradient increases southward as the zone of accommodation between the GB and WL narrows, consistent with previous studies (Kreemer et al., 2012).
The component of velocity parallel to the profiles (Figure S4 in Supporting Information S1) reveals smaller trends that vary with latitude, which reflect differences in the azimuth of the Sierra Nevada/Great Valley microplate (SNGV-Figure 1 and Figure S1 in Supporting Information S1) motion from south to north.The change in sign of the trend is related to the rotation of the SNGV microplate, which results in motion away from the Pacific plate in the south, and motion towards the Pacific plate in the north, consistent with SNGV counterclockwise rotation (Figure 1b) (e.g., Argus & Gordon, 1998;Bennett et al., 2003;Dixon et al., 2000).This rotation sets the western boundary condition and exerts a strong control on the strike and style of faulting east of the SNGV in the WL, which varies from south to north (Wesnousky, 2005;Wesnousky et al., 2012).

Strain Rates
Maps of the tensor strain rate components of shear and dilatation depict the geographic variation of strain accumulation which may eventually be released in earthquakes.In the fault slip rate analysis described below, we use a strain rate map to provide a regularization to the block modeling.Various kinds of parameterizations and regularizations for building strain rate maps have been tested and compared (e.g., Beavan & Haines, 2001;Hearn et al., 2010;Maurer & Materna, 2023;Shen et al., 2015;Spakman & Nyst, 2002;Tape et al., 2009).Our formulation separates gradients in the GPS velocity field on a sphere into components of tensor deformation and rotation using the parameterization of Savage et al. (2001).At each grid point we estimate the local strain rate tensor with a linear inversion from the gridded horizontal velocities as input data weighted with the velocity uncertainties.The results are insensitive to outlier velocity data and irregularities in station spacing because the velocities have been spatially filtered and gridded.Based on the result of an analysis of trade-off between high spatial resolution and low data misfit (Figure S5 in Supporting Information S1) we adopted a length scale of 8 km for the shear and dilatational component stain rate maps (Figure 3).We depict the shear as the difference between principal strain rates ( ė 1 ė 2 ) and dilatation as the sum of principal strain rates ( ė 1 + ė 2 ) (Figure 3).
The contiguous band of high strain rates east of the SNGV is the geodetic signal of the WL, its intensity decreases from south to the north.Its lower intensity in the Northern Walker Lane (NWL) is a function of the wider and shallower gradient in GPS velocity compared to the Southern Walker Lane (SWL) profiles (see B, C, D and E in Figures 1 and 2).The total velocity difference across the system decreases from ∼10 mm/yr in the SWL to ∼7 mm/yr in the NWL.Hence strain rates increase southward as a larger total velocity budget is accommodated across a narrower zone.
We discuss other features of the WL strain rate maps in parts: In the Southern WL (SWL) there are significantly higher strain rates in the Owens Valley (>80 × 10 9 /yr) adjacent to the eastern SNGV than there are along other fault systems of the eastern WL such as Death Valley, Furnace Creek, Fish Lake Valley, and Panamint Valley (see Figure S1 in Supporting Information S1 for map with fault names).This is attributable to the higher gradient in GPS velocities adjacent to the SNGV (see profile D of Figure 2).Strain rates across the Owens Valley are strongly transtensional, as seen in the relatively high shear strain rate (Figure 3) and positive dilatation rates (Figure 4).This is consistent with the strike of the Owens Valley (∼N20˚W) being in a releasing geometry with respect to azimuth of SNGV motion (∼N50˚W) (as noted by Unruh et al., 2003).The higher strain rates are more concentrated in the northern sections of the Owens Valley Fault and widen to the south to distribute deformation more evenly across the southern Sierra Nevada, Little Lake, Airport Lake, Panamint Valley, and southern section of the Death Valley Fault.The Death Valley/Fish Lake Valley fault system has lower strain rates (30-60 × 10 9 /yr) but are still elevated above the rates in the GB east of the WL.In the Central WL (CWL) the location and width of the zone of highest strain rates varies with latitude.Just north of the SWL, the shear steps left and concentrates, passing through the Long Valley Caldera near the SNGV (shown in detail in Hammond et al., 2019) and then widens again northward.Here the higher shear strain rates occupy all the area between the SNGV and northwest striking dextral fault systems east of Walker Lake such as Petrified Spring, Benton Spring, and Gumdrop Hills faults (Angster et al., 2020;Wesnousky, 2005).However, shear rates are highest at the east and west margins of the CWL, near the SNGV and east of the Wassuk Fault.However, dilatation rates (Figure 4a) show faster extension near the western margin adjacent to the SNGV and the Wassuk Fault.This is consistent with the partitioning between shear and extensional domains (Surpless, 2008).However, it may be only a matter of degree because the shear rates are lower west of the Wassuk but not zero, and the analyses of Bormann et al. (2016) suggests that some strike slip through this part of the CWL is required to explain the GPS data.At the left lateral faults in the Mina Deflection such as the Rattlesnake, Excelsior, Candelaria Hills faults the shear strain rates are lower (40-60 × 10 9 /yr) than those closer to the SNGV, but are still higher than those in the GB to the east, where the strain rates drop precipitously to below 10 × 10 9 /yr.
In the Northern WL (NWL) the strain rates separate and focus into two main arms, one in the westernmost NWL overlapping the Lake Tahoe Basin adjacent to the SNGV, and the other extending north-northeast to follow the Fairview, Dixie Valley and Pleasant Valley faults.This zone is the location of a sequence of earthquakes that released most of the seismic moment in the GB in the twentieth century, known as the Central Nevada Seismic Belt (CNSB).In the next section we discuss the transient nature of this active deformation feature.North of Lake Tahoe the NWL shear strain rates are lower (∼40 × 10 9 /yr) than to the south (>80 × 10 9 /yr) and confined to a narrow zone between the Pyramid Lake and Mohawk Valley Faults, narrowing even further north to lie between the Eagle Lake and Almanor Faults.Further north still a pocket of medium-high strain rates (∼50 × 10 9 /yr) lies at the Inks Creek, Battle Creek, and Bear Creek faults, a place that may represent the locus of active convergence between the SNGV and crustal blocks in the southern Cascadia forearc (Angster et al., 2020;Unruh & Humphrey, 2017).
Along the east edge of the WL eastward protrusions of elevated strain rates occur near locations of east-northeast striking left-lateral faults such as in the Carson Domain near the Olinghouse Fault, Mina Deflection, Gold Mountain east of the Death Valley Fault, and the Rock Valley Fault south of Yucca Mountain.The Clayton Valley and Lone Mountain faults also have rates slightly elevated compared to the rest of the GB east of the Death Valley fault, consistent with the findings of Lifton et al. (2013) who detected extension normal to strike near the level of the uncertainties with campaign GPS.The evolution of these faults may be linked to the nearby strike slip systems of Fish Lake Valley/Death Valley (Oldow et al., 1994), suggesting that perhaps all these protrusions are similar cases with similar mechanical origin.Significant strain rates protruding east of the WL, elevated above background GB rates have been seen in other recent strain rate models as well (e.g., Zeng, 2022b), and are especially apparent in maps made with combinations of geodetic and geologic data (Kreemer & Young, 2022).These signals may indicate that complexities exist in the deformation field in the GB east of the WL but are currently near the resolution of the geodetic data, which may be clarified with better geodetic coverage east of the WL and/or longer time series.
Inside the SNGV strain rates are very low (mostly <10 × 10 9 /yr), consistent with previous geodetic studies (Argus & Gordon, 1998;Bennett et al., 2003;Dixon et al., 2000;Kreemer et al., 2012), lower levels of seismicity and sparsity of faults in the compilations (Figure 1c and Figure S1 in Supporting Information S1).The boundary between the low strain rates in the SNGV and the high strain rates in the WL is very near the faults bounding the SNGV east edge.The most significant exception to the SNGV rigidity is in the southernmost part of the microplate west of the SWL and Owens Valley, where the Kern Canyon, Lake Isabella and western extension of the White Wolf faults extend north of the Garlock fault into the High Sierra.Seismicity in this area has been used to suggest that the southernmost Sierra is subject to heterogenous extension and crustal thinning associated with foundering of lower lithosphere (Unruh et al., 2014).
The dilatational strain rate field indicates rates of area change that are generally lower than shear strain rates (Figures 4a and 4b) but are positive on average consistent with the WL and GB being a transtensional part of the PA/NA plate boundary system.Owing to smaller signal-to-noise ratio dilatation rate is harder to resolve than shear rate in the WL, but some clear patterns emerge.Rates are faster in the WL that in the GB, and the stripe of faster extension adjacent to the Sierra Nevada coincides with a zone of extension-dominated transtension noted by Oldow (2003).Strong positive dilatation is observed near the Long Valley Caldera active magmatic system that is undergoes episodic inflation (Montgomery-Brown et al., 2015).This inflation affects a volume of crust with active strain penetrating westward into the SNGV and eastward into the WL, influencing seismicity which is transient deformation possibly not representative of long-term deformation (Hammond et al., 2019).

Correction for Viscoelastic Postseismic Relaxation
Some signals in the geodetic data are not representative of the long-term rate and pattern of deformation and must be removed from the velocity field before estimating time-invariant slip rates.In the WL and western GB signals from post-earthquake relaxation related to the late nineteenth and twentieth century M 6.9-7.5 earthquakes in the CNSB (Bell et al., 2004;Wallace, 1984) are apparent in the geodetic data.It has been shown in several studies that these earthquakes initiated uplift and dilatation anomalies in central Nevada, which stand out from the slow background rates in the rest of the Basin and Range (Gourmelen & Amelung, 2005).We use a model of the viscoelastic relaxation process for the CNSB (Hammond et al., 2009) using the theory and software of Pollitz (1997) to subtract the transient anomaly from the measured velocities.The model accounts for the effects of nine earthquakes that occurred east of the SNGV and within our area of study including the 1872 Owens Valley M7.4, 1915Pleasant Valley M7.3, 1954Dixie Valley M6.9, 1954 Fairview Peak 7.0.The strain rates predicted from the relaxation model have concentrations of both shear and dilatation owing to the obliquity of the ruptures and uniaxial component of extension at the CNSB.The predictions illustrated in Figure S6 in Supporting Information S1 were interpolated to the same GPS stations used in this analysis, subject to the same GPS Imaging process before tensor strain rates are calculated, and mapped with the same bounds as Figures 3 and 4.There are uncertainties in the relaxation model associated with limits of our knowledge of the slip during these past earthquakes as well as in the structure and rheology of the WL crust and mantle.However, the model is based on seismic and geologic data and estimates crust and mantle viscosities similar to those found in other studies in the western United States.Details of model construction are provided in Hammond et al. (2009) which includes description of the rheological assumptions that were made, how it is derived from earthquake parameters.Another indication that this model is reasonable is that the strain rate field corrected for CNSB postseismic relaxation removes the long finger-shaped anomaly of high shear and dilatational strain rate that branches north-northeast from the main band of high strain rates northeast of the CWL (Figures 4a and 4b).In the corrected strain rate maps the eastern boundary of the higher strain rates in the WL becomes much straighter and more parallel to its western boundary (Figure 3b).
While the model removes the single most obvious north-northeast trending anomaly, we notice the existence of several other lower intensity anomalies that extend in a similar direction to the CSNB anomaly, but are not corrected by the relaxation model.One is near Pyramid Lake, NV and the other near Eagle Lake, CA, west of the CNSB.Whether these reveal the unmodeled effects of earthquakes occurring further in the past is unknown.However, it is likely that not all events occurring in the last 1000 years on WL faults are documented.Paleoseismic data show, however, that an event with ∼3.3 m of slip occurred on the Incline Village Fault in the northern Lake Tahoe Basin ∼500 ± 150 years before present (Seitz et al., 2016).The amount of offset and length of the fault suggests a magnitude of ∼7.1.That would have been enough seismic moment to generate a postseismic viscoelastic signal that is detectible with GPS, as in the CNSB.Whether that signal persists to the present day and is related to the observed strain rate anomalies is speculative.There are other low amplitude strain rate anomalies in both the shear and dilatation fields that are noticeable east of the WL, where background strain rates are lower.However, some of these may be related to artifacts in the imaging technique.
We do not attempt to correct the GPS velocity field for the magmatic inflation at LVC because volcanic deformation is not a cyclic effect as is the case for faults.Inflation-related deformation may be cumulative and influence nearby faults so should be modeled as part of the strain accumulation field, and related slip rates.

Justification for the Approach
Block modeling assumes that the crust is divided into contiguous moving elastic volumes that drive slip on the faults that bound them.The method accounts for the fact that the data are collected during the interseismic period when fault systems are locked from the surface to the bottom of the seismogenic upper crust but slip at the block motion rate at greater depth.The parameterization enforces kinematic consistency between fault slip rates and crustal block translation and rotation.Several implementations have been developed, varying in parameterization, complexity, and use of regularization (e.g., Bennett et al., 1996;Evans et al., 2015;Loveless & Meade, 2010;Matsu'ura et al., 1986;McCaffrey, 2002;Meade & Hager, 2005;Savage & Burford, 1973).
In the WL a system-scale block modeling approach is particularly needed because there are many faults that together form network which creates a complex distribution of seismic hazard.Block modeling represents the complexity in the network of faults with a framework that allows integration of several kinds of data to understand the system.Block modeling is, however, challenging for a couple of reasons.First, vertical-axis rotations of crustal blocks are observed in the WL and these rotations are intimately linked to the fault slip.While paleomagnetic data constrain rotation since 10-13 Ma for some blocks (Carlson et al., 2013;Cashman & Fontaine, 2000;Petronis et al., 2009), their present-day rates of rotation are not well known, and data is not available for all WL areas.Given the spacing between faults (5-30 km) and locking depths (∼15 km) the signals of elastic strain accumulation and rotation across the blocks spatially overlap.This leads to the solution for block rotations and slip rates to be under-determined, with trade-offs between parameters.Regulation of the problem is required, and our approach to this is explained below.Second, the large number of faults and complexity of the fault system (Figure 1c) make manual building of traditional block models cumbersome.Block models generally require completely connected boundaries that define independent polygonal domains.However, it is not always possible to define the boundaries objectively because mapped fault traces are discontinuous, may be based on incomplete data sets, and must be drawn according to subjective decisions from the analyst.We address this difficulty in the next section.

Automated and Iterative Block Model Generation
Here we describe the procedure for generating block models from a database of faults.To begin we require that each fault in the database be represented by an ordered sequence of coordinates that define the surface trace and the fault dip.We use the NSHM database of western US faults (Hatem et al., 2022a) which provides these parameters.There are 373 faults in the database that touch the study domain for which we estimate a slip rate (Figure 1).While not every known fault in the GB is represented in this database, it includes the structures that are best studied and have demonstrated Quaternary activity.
Prior to constructing each block model the fault traces are further simplified beyond the already simplified traces available in the NSHM database (Hatem et al., 2022c).We select a subdomain consisting of a latitude/longitude box of 2°× 2°and truncate the fault traces to include only those segments inside the box (Figure 5).We also down-sample the number of nodes used to define fault traces so that each segment has length no smaller than ∼ one fortieth of the dimension of the subdomain (0.05°).The spacing of segment nodes is controlled because they define seed locations for the next step which is the generation of an initial set of blocks.Nodes on the rectangular boundaries of each model are needed so that blocks can extend to the edge of the domain.We place one at each corner and 3 more along each side so that the initial configuration has blocks that are small enough to represent fault geometries.The result is a set of nodes on the faults and boundaries that define a Delaunay triangulation where every fault segment is the side of one of the triangles.Associating each fault segment with a triangle side makes a valid block model that honors the fault network (Figure 5b).However, this model is primitive because the presence of many small and narrow blocks makes it poorly parameterized for deformation modeling.The model suffers from an excessive number of free parameters (3 Euler rotation vectors for each block) and the block rotations are difficult to constrain with GPS velocity gradients because block dimensions may be small in one or more directions.
We reduce the number of blocks in the model by joining selected abutting block pairs into single blocks and repeating the process iteratively.Since each block is represented by an ordered sequence of node numbers, two blocks are abutting if they share a subset of nodes.Removing the boundary between two blocks involves deleting them and then adding a new one composed of the polygonal union of the original two.At each iteration we prioritize which block pairs to combine by developing a ranking system, where block bounding segments are scored according to their geometric properties.Segments with low scores are selected for elimination, and the two adjacent blocks are joined into one.Block boundary segments that align with faults are given score of infinity so they cannot be lost through joining of blocks.Segments that are aligned with the domain boundary are also given score of infinity.Segments that otherwise touch faults are given a score equal to the number of fault segments touched, which preserves some detail in the model near faults, but also allows for simplification so that individual segments along a fault tend to have similar slip rates.Blocks are also targeted for combination if they have very small angles at vertices (20°or less), small area (100 km 2 or less), large shape parameter defined as the perimeter divided by the square root of the area (10 or more, where a circle is ∼3.5), or large interior angles (190°or more) which prevents blocks with significant concavity.Also blocks having more than one contiguous chain of boundary nodes ("donut blocks") are forbidden by not allowing any block merging that results in this condition.An example of a block number reduction operation sequence is shown in Figure 5.The systematic joining of many small blocks into fewer larger blocks (e.g., from Figures 5b and 5c) results in a model that is better conditioned for the inversion procedure because there are fewer model parameters to be constrained by the available geodetic data.The procedure ensures that the important parameters, that is, slip on real geologic faults in the database, are preserved.After a round of block number reductions, the process is repeated until the number of blocks is small enough (less than 50 blocks), the number of blocks does not change, or the algorithm fails by generating a block model that is invalid.The most common reason a model becomes invalid is when the combined block is not correctly parameterized with an ordered sequence of nodes shared with an adjacent block.This condition can arise (rarely) owing to errors during the combination operation.Because of the large number of models generated, loss of a few owing to block combination errors is acceptable.
Our strategy for increasing robustness in the slip rate estimates is to generate many models to reduce dependence on any single model's representation of the fault network.While the generation of each model is deterministic (the same fault inputs will always result in the same block model geometries), small changes in the location of the bounding box with respect to the faults result in a different initial node configuration and Delaunay triangulation.All models have block boundaries that coincide with the input faults, but the other off-fault block boundaries traverse the area between the faults in different places.We generate a grid of overlapping rectangular model subdomains with a spacing of 0.25˚× 0.25°.At each iteration the subdomain boundary is shifted by 0.25°which changes the relative position between the rectangular boundary and the nodes along faults inside, resulting in a different triangulation.The grid extends beyond the area of interest (Figure 1), but if the intersection between a sub-domain and the full model domain (Figure 1) is <10% of the sub-domain area then the model is not used.A total of 1,240 block models are built using the generative procedure described above and 1,231 of them (99%) were valid and could be used to estimate slip rates from the GPS velocity field.

Solving for Slip Rates
For each block model we estimate the slip rates on each fault segment using the method of Hammond et al. (2011), which assumes the geodetic data are collected during the interseismic time when the faults are locked at the surface but the block motions continue.This method is like others (e.g., McCaffrey, 2002;Meade & Hager, 2005) but uses a damped linear least squares inversion to solve simultaneously for a set of block rotations and fault slip rates.The motion of each block j is parameterized with an 3 × 1 Euler rotation vector ( ω with parameters representing strike-slip (a k ) and dip-slip (b k ) rates for fault k.These unknown parameters are related to the horizontal vector GPS velocity data v → GPS,i at station i through the equation: where the first term on the right side is the contribution from block rotation at the station location r i , and the second term accounts for elastic strain accumulation owing to shallow locking of the faults.The functions G ss,ki for strike slip and G n,ki for normal/thrust slip are the Green's functions that enforce back-slip (Savage, 1983) and are calculated using functions from Okada (1985), which are fixed by the known location of the faults and stations.Multiple faults may contribute to elastic strain accumulation at a single station, up to a maximum of L faults, which we set in this analysis to 6 since greater values do not change the result significantly.
Regularization is implemented by minimization of both fault slip and the vertical axis spin rates.These conditions are expressed as: ) which are equations added to the system for all k and j.Adding the constraints in Equation 2to the inversion problem enforces minimization of the slip and rotation rates while also fitting the data.The weight of the damping condition is controlled by the regularization values that are a function of background strain rate and style.This is needed because of the orders of magnitude variation in strain rate in different locations in the WL makes spatially constant factors ineffective in some areas.Strike slip and normal/thrust components have damping weighted separately so that slip rates are consistent with the tensor strain rate style.The regularization method is discussed in detail in the Supplemental Materials.In accordance with previous studies of the WL seismogenic depth (e.g., Ruhl et al., 2020;Zuza & Cao, 2020), in each model we assume 15 km locking depth throughout the entire area.
No parameters for long term strain rates within individual blocks are included.
The median slip rates for each fault from the set of all block models is shown in Figure 6 and are provided in Supplemental Table S2.Taking the median value reduces the impact of outliers that can be seen in solutions of some individual block models (e.g., Figure 5b).Each fault has multiple segments with some variability along strike so we consider each segment as an individual sample of the slip rate which contributes to the distribution.Faults with a strike slip component have a median number of 454 individual segments contributing to the slip rate estimate (over all models and segments), but as few as 12 and as many as 3036.Long faults tend to have more individual segments, and faults near the boundary have fewer since not as many block model sub-domains cover them.There are fewer dip slip rates estimated because we did not estimate them for faults which are deemed to be vertical strike slip faults in the NSHM database.Some block models produce a slip rate that is inconsistent with the geologically determined sense of slip for the fault.We truncate the distributions of slip rates in all cases where the slip rate has a sign that is in direct conflict with the geological characterization of the fault rake and determine the median rates from the remaining distribution.East of the WL, in the Basin and Range, most of the faults were designated as normal faults in the NSHM database, so had no strike slip rate estimated.However, the normal rates were estimated and had diversity that is not possible to appreciate given the 1 to 1 mm/yr color scale limits in Figure 6.To see the diversity of normal slip rates in the Basin and Range we provide another version of Figure 6b in the Supplementary Materials with color scale to 0.1 to 0.1 mm/yr (Figure S7 in Supporting Information S1).

Slip Rates
For discussion we use the convention where dextral slip has negative and sinistral has positive sign, normal slip has negative and thrust positive rates.For a given fault we report the median slip rate for all segments and all block models that contribute an estimate.We estimate the uncertainty with the median absolute deviation (MAD) of the slip rate estimates, multiplied by 1.4826 which makes the value identical to the standard deviation if the distribution of estimates is Gaussian (Wilcox, 2005).In the case of the automatically generated block models there are occasionally instances when the model returns an outlier slip rate that is not representative of the body of the distribution.Thus, using the median and MAD is advisable because they are insensitive to outliers and more representative of the body of the estimates.The resulting rates and uncertainties account for variability in fault strike and in the geometric depiction of the local block kinematics.
We do not divide by the square root of the number of model estimates as is done when estimating the formal uncertainty in a mean because the large number of models generated (sometimes thousands-Table S2) would make the uncertainties unrealistically small.However, faults that have high degree of strike variability or higher slip rates will tend to have larger uncertainties not because of greater error in the models, but because the geometric variability along the fault maps to variability in slip rate, or their geometries are difficult for the block generator to represent consistently.The resulting set of fault slip rates reveals the distribution of active deformation into domains in the major sub-provinces of the WL (Figure 6).
Strike slip rates are generally higher in the SWL, CWL and NWL compared to the rates in the Basin and Range east of the WL.For example, in the SWL the long northwest striking dextral faults accommodate the largest amount of deformation.The fastest dextral rate is 4.2 ± 1.4 mm/yr (dextral) for the Hunter Mountain/Saline Valley Fault, followed next by 2.9 ± 1.0 mm/yr (dextral) for the Death Valley Black Mountains segment, 1.5 ± 0.6 mm/yr for the Owens Valley Fault, and 0.8 ± 0.6 for the Panamint Valley fault.These faults account for most of the velocity budget of across the SWL (Figure 2).While there is a strong contrast between high slip rates in the SWL and low rates in the Basin and Range exclusive of the WL, some SWL slip rates are low, even though they reside in areas with higher strain rates.For example, the extension rates of the Deep Springs Valley ( 0.4 ± 0.1 mm/yr), Dry Mountain ( 0.2 ± 0.2 mm/yr), Tin Mountain ( 0.1 ± 0.1), Towne Pass ( 0.6 ± 0.3) faults are systems that accommodate northwest-southeast oriented extension amid the ranges that lie between the Death Valley, Panamint Valley, and White Mountain fault systems.
In the CWL there is evidence of partitioning, where faster strike slip rates (faster than 1 mm/yr) step right at the Mina Deflection from Fish Lake Valley and White Mountain into the Gumdrop, Benton Springs and Bettles Well/ Petrified Springs fault systems east of Walker Lake.This is consistent with geological observations of fault rakes (Angster et al., 2019;Pierce et al., 2021;Surpless, 2008;Wesnousky, 2005) and in our model is in part a consequence of imposing constraints on rake from the geologic database (Hatem et al., 2022a), which suppresses strike slip on faults with rake = 90.The relatively high dextral slip rates stay to the east edge of the WL until further north in the NWL, north of the latitude of Lake Tahoe, they become distributed across the width of the WL between the SNGV and Pyramid Lake (Figure 6).
In the NWL the higher dextral rates associated with the WL extend at least as far north as the northern end of the SNGV microplate, including the Lake Almanor, Eagle Lake and Likely fault systems (faults names are labeled on Figure S1 in Supporting Information S1).We see thrust sense on faults immediately north of the SNGV microplate, where it interacts with the Klamath Mountains and Oregon Coastal microplate at the Red Bluff, Bear Creek and Inks Creek faults (Figure S1 in Supporting Information S1), consistent with expectation based on other geodetic, geologic, and seismic studies (Angster et al., 2020;Unruh & Humphrey, 2017;Williams et al., 2006).The rate on the Hat Creek Fault to the north is slow and mostly normal sense ( 0.1 ± 0.1 mm/yr).Here, normal slip was not permitted on these faults since in the geologic database their dips were 90°.
The model images sinistral slip along the east-northeast striking faults that cross the WL and reach to its eastern edge such as the faults of the Mina Deflection, Garlock and Olinghouse faults, and Carson lineament.Sinistral slip on these structures is consistent with median model clockwise vertical axis rotations that are discussed below, with expectation based on models of blocks rotating in northwest directed shear field (Bormann et al., 2016;de Lano et al., 2019;McKenzie & Jackson, 1983;Platt & Becker, 2013), and with observations of paleomagnetic clockwise rotation (Carlson et al., 2013;Cashman & Fontaine, 2000;Petronis et al., 2009).Sinistral slip is seen at the southern end of the SNGV at the White Wolf (4.2 ± 1.2) and Garlock Faults (0.3-1.2 mm/yr increasing to the west).The Garlock fault slip rates are slower than those estimated in most geologic studies, 2.7-5.3 mm/yr or higher, see Hatem and Dolan (2018) and references therein for a summary.We include only its central and easternmost segments here, which are slower than its western section, in both this study and in geologic observations (McGill et al., 2009).
Normal slip rates in the WL are also generally higher than normal slip rates in the Basin and Range to the east, which are all closer to zero than 0.2 mm/yr.While the WL strain rate field is dominated by shear, it also has enough positive dilatation and releasing bends in the fault systems to drive normal slip.Some long faults have relatively fast normal slip rates with uncertainties 0.5 mm/yr or less (e.g., northern Kern Canyon 1.4 ± 0.5 mm/ yr).The very highest normal slip rates are on shorter faults that tend to have high uncertainties (e.g., Hartley Springs at 1.8 ± 1.5 and Airport Lake 1.3 ± 1.0 mm/yr).These are both in locations with multiple intersecting or near-overlapping faults, so the block construction algorithm may be drawing boundaries around these faults in a greater variety of ways.Extension occurs on most faults north of the Garlock if they do not dip 90°in the fault database, on which normal slip was not estimated (gray faults in Figure 6b).In the CWL normal slip is distributed between Lake Tahoe and the Wassuk Fault near Walker Lake.In the area northwest of Lake Tahoe, between the Mohawk Valley Fault and Pyramid Lake has positive dilatation, but the faults in the NSHM database are all vertical strike slip faults, so normal slip was not estimated.
Normal slip on the Lone Mountain (0.6 ± 0.4 mm/yr), Clayton Valley (0.5 ± 0.4 mm/yr) and Emigrant Peak (0.4 ± 0.3 mm/yr) faults, south of the Mina Deflection and east of the WL show extension across the Silver Peak/ Lone Mountain complex and appear consistent with displacement-transfer style faulting (Oldow, 1992;Oldow et al., 1994) and geological fault slip rates to within uncertainty (Foy et al., 2012;Lifton et al., 2015).
In the Basin and Range east of the WL the median normal slip rate is 0.05 mm/yr and robust measure of their standard deviation is 0.05 mm/yr.This suggests that the geodetically measured deformation observed across the eastern Nevada Basin and Range (Hammond et al., 2014) is discernible in the slip rates.It is also consistent with the strain release rate observed in paleoseismic trenches on the active normal fault systems between 38.5°and 40.0°that cumulatively add to between 0.8 and 1.0 mm/yr of extension across 450 km (Koehler & Wesnousky, 2011).If ranges are separated by ∼30 km then 450 km equates to 15 ranges or 0.07 mm/yr extension per range, similar to our median geodetic normal slip rate.
While Hammond et al. (2014) found a shear sense of strain rate in the geodetic data, there is no systematic strike slip component in the slip rates estimated here or in the trenches of Koehler and Wesnousky (2011).While the signal of extension is clear in the normal slip rates, there is a possibility that our regularization based on scaling by our GPS strain rates model could be suppressing the signal of strike slip since the strain rate in eastern Nevada is very low (<5 × 10 9 /yr).Another possibility is that the shear observed in geodesy is a transient deformation associated with the sum of many past and possibly distant western US earthquakes.This is predicted by the forward model of Young et al. (2023) at a level below 2 × 10 9 /yr.Removing their predicted postseismic strain from the observed strain rate makes it more uniaxial with an east west extension direction similar to the geologic extension direction.Whether transients are affecting the GB in this way deserves more direct study of the potentially far-reaching impact of post-earthquake viscoelastic transients.

Slip Partitioning
We provide an example that shows the effectiveness of the method and regularization in partitioning transtensional strain appropriately between closely adjacent strike slip and normal faults.The Owens Valley fault is a vertical strike slip fault and is closely adjacent to (3-15 km) the Independence Fault which is an east dipping normal fault that bounds the Sierra Nevada range front (Figure 7).Their proximity means that the motion of the area between them must be constrained in a robust way to accurately estimate both slip rates in a block model.Adding to the challenge is the asymmetry in the strength of GPS data coverage, with strong network coverage in the SWL to the east and weak coverage in the Sierra Nevada wilderness areas to the west.Our method finds a strike slip rate for Owens Valley of 1.5 ± 0.6 mm/yr (dextral) and for the Independence Fault 0.1 ± 0.1 mm/yr for its strike slip rate and 0.5 ± 0.30 mm/yr normal slip rate.These rates are in accordance with geologic observations in terms of rates and style (Bacon & Pezzopane, 2007;Beanland & Clark, 1994;Haddon et al., 2016;Jayko, 2009;Lee et al., 2001).While the normal component of slip on the Independence Fault is only significant at the 1-sigma level, the higher level of uncertainty is to be expected here because of the lack of GPS constraint west of the fault and the variability of its strike.This example suggests that the combined power of the features of the method results in robust estimates of fault slip rates that abide by geologic constraints, provide realistic uncertainties, without explicitly constraining the rates to geologically determined rates a priori.

Vertical Axis Block Rotations
Rigid block motion on the surface of a sphere can be decomposed into two orthogonal components.The first is translation of the block associated with an Euler rotation vector 90°distant from the block centroid, the second is vertical axis spin which is associated with an Euler rotation vector parallel to the direction pointing to the block centroid from Earth center.For each block we separate its estimated Euler rotation vector into these two components and then using the vertical axis spin component estimate at each pixel on the map the median spin rate from all the models (Figure 8).
The resulting image shows that in the WL spin is mostly clockwise with values less than 2°/My.
Pockets of faster spin rates occur near the locations of east-northeast sinistral faults, including the Mina Deflection, Olinghouse Fault, Carson lineament, Garlock Fault and faults south of it.The association between sinistral faulting and clockwise rotation is consistent with the paleomagnetic data and models of rotations as noted in the previous section.There is also clockwise rotation in the CWL west of the Wassuk Range Fault and east of Smith Valley of ∼1°/My consistent with other recent CWL block models (Bormann et al., 2016).
Even greater degrees of vertical axis spin rate are seen north of the SNGV.However, in this area the geographic density of faults in the NSHM sources database is very low.For example, there is an usually large gap between the Pondosa proxy fault and the next fault to the west which is the Trinidad Fault near the California coast 200 km west (off of map in Figure S1 in Supporting Information S1).Thus, clockwise spin is needed to accommodate GPS velocity gradients that are part of the regional clockwise rotation pattern between the SNGV and Oregon Coast microplate (Hammond & Thatcher, 2005;Unruh & Humphrey, 2017;Williams et al., 2006).Also, there is an unusually large gap of 74 km between the active Surprise/Warner Valley faults and Steens Mountain fault systems (Personius et al., 2007(Personius et al., , 2009) ) in which there is no fault in the sources database, though dozens of shorter faults are present in the USGS QFFD (USGS and AGS, 2011;USGS and CGS, 2011;USGS and NBMG, 2011;USGS and UGS, 2011).Here there are moderate clockwise rotation rates in our model ( 0.5 to 1.0°/Myr), possibly because of fewer faults are present in fault database so the GPS velocity gradients tend to be accommodated through block rotation rather than fault slip.If faults missing from the database are included in a future release, the rates of block vertical axis spin needed to explain the GPS data may be reduced.Crossing southern Nevada is a zone of positive (counter-clockwise) spin that follows a band of seismicity that extends from southwest Utah to the SWL.Known as the Southern Nevada transverse zone (Slemmons et al., 1965) it has been characterized with GPS measurements as a zone of sinistral deformation transfer between the Wasatch Fault system and SWL called the Pahranaghat Shear Zone (Kreemer et al., 2010).However, this band is in a zone where there is a low density of fault segments in the database, and so high rotation may be imaged.Also, this zone has lower GPS station density (Figure 1b), making it more difficult to resolve the zone's location precisely.It does extend all the way to the Death Valley fault, and crosses it to the Hunter Mountain/Saline Vallely fault, encompassing Tin Mountain and Dry Mountain faults.

Model Uncertainties and Data Misfit
Misfit of the models to the data is obtained by finding for each GPS station the median predicted east and north velocity from the subset of block models that use the station.The residuals are the difference between the original GPS velocities (not the gridded or filtered velocities) and the predictions (Figure 9).The figure annotation shows the MAD times 1.4826 of residuals is significantly smaller than the standard deviation indicating that outliers affect the standard deviations substantially (by about a factor of 3 over the robust estimate made with the MAD).The histograms of residuals normalized by their uncertainties have MAD*1.4826 of ∼2 indicating that the data are fit at a level about twice the uncertainties in the velocities, with a robust estimate of RMS of ∼0.7 mm/yr.Our misfit is a bit lower than those in the NSHM western US models (Johnson et al., 2024) in part because we use the robust estimate which reduces impact of outliers and we focus on a subset of the western US with slower deformation.
The histograms show near zero mean for the east and north residuals, with a slight tendency for positive mean north residual at the level of ∼0.1 mm/yr.To address the potential for systematic misfit we use a definition of systematicity that for each station takes the mean dot product of neighboring station velocities (defined by Johnson et al., 2024): We use a radius of 30 km in the vicinity of each station i to select the n nearest stations.When velocities have similar azimuths their dot products have higher magnitudes and S i increases, while when the vectors have random azimuths the signed dot product values tend to drive S i toward zero.This measure varies between 0 for very randomly oriented vectors to 1 when they all have the same azimuth.The absolute value of S i for each station is plotted in Figure 10a, and reveals domains with significant systematicity, similar to the models that comprise the NSHM deformation model suite (Johnson et al., 2024).However, this measure does not consider the uncertainties in the velocities so may give the visual impression that large areas have systematic residuals, even if the residual velocity magnitudes are below the uncertainties in the data.
To address this we define an alternative version of the formula that normalizes the residual vectors by their uncertainties instead of the residual norms: When we look at Figure 10b which is based on Equation 4we see that many of the residuals that were highly systematic in Figure 10a have residuals that are well below the uncertainties in the velocity data.This is expected because these small residuals will not influence the least squares inversions for block rotation and slip rate parameters as much as other more significant residuals.There are, however, still a few significant systematic residuals in the vicinity of Long Valley Caldera where magmatic inflation drives radial signals that are not well modeled on the CWL fault systems.Also, in the NWL near the epicenters of the 1954 Fairview Peak and Dixie Valley earthquakes, where the CNSB postseismic relaxation has been modeled and removed (see Section 3.2), we see an indication of systematic misfit with Sʹ ∼3-4.This suggests that the signal from CSNB postseismic relaxation may not be completely removed and some signal that is not well modeled by CNSB faults remains.The third area with high Sʹ values is east of the dextral faults in the eastern CWL near the Toiyabe Range.This anomaly is near the source of the 1932 M7.1 Cedar Mountain earthquake, which is another events included in the CNSB transient postseismic relaxation model and so may also be indicating that the transient signal is not entirely removed.

Comparisons Between Geologic and Geodetic Slip Rates
Geologic and geodetic slip rates each aim to measure the same thing, that is, the relative rate of motion of blocks of Earth's crust on either side of a fault.However, they use very different means to measure displacement over very different time scales.Geologic data constrain relative motion across faults and in landscapes over the long periods of time needed to accumulate measurable displacements of earth materials.Geodetic data, because of their high precision can constrain displacements over a few years to decades (Bennett, 2007).Similarity between geologic and geodetic slip rates indicates that the different methods are estimating the same potential for slip at different times in the earthquake cycle, which increases confidence in the results.Because strain rates in the WL are slower than along the main plate boundary fault zones in western California, our slip rates occupy a much slower and narrower range than in other recent comparisons.It is therefore more challenging to establish a correlation between geodetic and geologic rates in the WL.It is important to note that unlike in some other studies (e.g., Pollitz, 2022;Shen & Bird, 2022;Zeng, 2022aZeng, , 2022b) ) we do not use geologic slip rates as prior constraints or bounds on the solution for slip.Thus, our geodetic rates have a higher degree of independence from the geologic rates and are corroborative of them when they agree.
We compare our geodetic slip rates to geologic rates from a recent compilation developed to support the NSHM (Hatem et al., 2022b, which we also provide in Table S2).We exclude faults south of 34.8°latitude and the Garlock fault because its geologic rate in the NSHM database (up to 11 mm/yr) is much greater than the range of all other faults in this study and our estimate (1.2 mm/yr in the central section).The result (Figure 11a) shows a degree of agreement but many slip rates are less than 1 mm/yr so we also plot them on a log 10 axes (Figure 11b).The log scale plot reveals a general trend of agreement between the geodetic and geologic slip rates.The correlation between slip rates is similar on linear (r = 0.37) and log 10 scales (r = 0.44).This is a weaker correlation compared to rates across the entire Western US, which are between 0.41 and 0.88 depending on the contributing modeler (Johnson et al., 2024).However, a lower correlation is expected given the far narrower range of slip rates in the WL (all are less than 5 mm/yr), and because we do not impose a constraint in our modeling that our slip rates should be near the geologic rates.Nonetheless our conclusions are similar to what was observed for the entire western US: that the faults with the lowest geologic rates (<0.05 mm/yr) have higher geodetic rates (Figure 11b).If we reckon that slip rates agree when the geodetic rates are within the geologic minimum/maximum bounds or when the geologic rates are within 2 times the geodetic uncertainties (to be within their 95% confidence interval), then 70 out of 362 slip rates (19.3%)disagree and the rest (80.7%) agree.
We also compare our geodetic rates to the median geodetic rates from the western US geodetic modelers (Figures 11c and 11d).The correlations are slightly higher than the correlation with the geologic rates, with r = 0.40 for linear scale and r = 0.42 for log 10 scale.Because we solve for the same parameters from the same data, we may expect some agreement with their results.However, many of the NSHM modelers optimized their analyses for (a) a much larger and more tectonically diverse area with slip rates that vary by several orders of magnitude, (b) used a variety of analytical techniques, and (c) used geologic rates as prior constraints.The similarity of correlation shows that our rates agree with other geodetic slip rates about as well as they agree with geological slip rates.This may possibly be because the NSHM modelers relied heavily on geologic rates.
In the GB east of the WL, many of the faults have low slip rates, both in our model and in the NSHM database.In this area many geologic rates in the NSHM database have slip rates exactly 0.1 mm/yr because they had no preferred rate, only a minimum of 0.0 mm/yr and maximum of 0.2 mm/yr.In these cases, we took the geologic rate to be the mean of the minimum and maximum (Figure 11b).Our geodetic rates for these faults vary between <0.01 mm/yr to over 1 mm/yr.The median difference between geologic and geodetic rates in this area is 0.03 mm/yr which is close enough to zero to indicate that there is no bias for geologic or geodetic rates being greater.The MAD of the differences is 0.10 mm/yr.The absolute level of disagreement is very small mostly because there are so many faults in the Basin and Range with very low slip rates, and any slip rate estimate consistent with the low background strain rate would be similar at the level of 0.1 mm/yr.It is an advantage of our robust multi-block method that it returns the geologically plausible slip rates for large areas of the GB with many low slip rate faults.
Lastly, in Figure 11e we compare the distributions of NSHM geologic and NSHM geodetic slip rates to our geodetic rates.The faults in each set are sorted from smallest to largest slip rate value and plotted in order.While it is difficult to show that our slip rates have greater accuracy compared to the true slip rates, the distribution shows that they have a more natural diversity in slip rate estimates especially below, and without clustering at 0.1 mm/yr.The similarity in distribution between the NSHM geodetic and NSHM geologic slip rates is likely a symptom of the dependence of NSHM geodetic rates on NSHM geologic rates.This may lead to a systematic overestimation of integrated moment in the Basin and Range if the 0.1 mm/yr value is used for all faults.

Off-Fault Deformation
Wesnousky et al. (2012) showed that there are paths that can be walked from the Sierra Nevada crest across the CWL that do not cross a mapped fault, suggesting that in some locations the relative motion between the SNGV and GB is accommodated without faulting.Moreover, the geodetic modelers for the western US NSHM explicitly quantified deformation in the GPS data that was not mapped onto faults and derived off-fault moment rates between 31% and 58% across the Western US.Off-fault deformation may represent strain in the deformation budget that occurs near, but not on the main fault strands (McGill et al., 2015;Oskin et al., 2007), or possibility entire crustal blocks undergo non-brittle deformation by developing folds or oroclinal flexures that accommodate strain (e.g., Faulds & Henry, 2008), or some other process.
We estimate on-fault deformation as moment from slip on the active faults in the NSHM sources database.We sum the moments by assuming a constant seismogenic thickness of H = 15 km, a shear modulus of G = 30 GPa, and use fault areas of the locked segments A = LW, where W = H/sin(dip) and L is fault length.Summing over all faults segments the quantity GHAs, where s is the slip rate gives a total on-fault moment of 4.3 × 10 18 N-m in our model domain.We estimate off-fault deformation by summing the same terms over all the boundaries between blocks that are not in the NSHM database.The block boundaries are not guaranteed to occur in the same place in every iteration or to have the same sign of slip.Thus, the off-fault moment is distributed within zones between the NSHM model faults (Figure 12).The total off-fault moment rate in the domain is 9.1 × 10 18 N-m, giving an offfault deformation proportion of 68%.
This proportion of off-fault deformation is higher than that found in western US models for the NSHM (Evans, 2022;Pollitz, 2022;Shen & Bird, 2022;Zeng, 2022a).There are two reasons why the proportion may be particularly high in the WL.First, we include explicit parameterization for off-fault deformation to occur, that is, on block boundaries that are not faults, and this allows the deformation to be detected and quantified.While the method allows more off-fault deformation to emerge, there is a limit to how well its location is imaged, since it is forced to lie where block boundaries not in the fault database are created by the automated block generation procedure.Second, in the WL and GB there are likely more faults undiscovered and/or not included in the NSHM database, even though they are active and contribute to the accommodation of far field budgets.Many faults in the USGS QFFD (USGS and AGS, 2011; USGS and CGS, 2011; USGS and NBMG, 2011; USGS and UGS, 2011) are not present in the NSHM database (Figure 12).Thus, the degree of completeness of the fault database may be lower in the WL and GB compared to other systems in the western US, for example, the San Andreas.A recent example is the 2020 M6.5 Monte Cristo Range earthquake in Nevada which ruptured the surface on unmapped segments of the Candelaria Fault (Koehler et al., 2021).
The map of off-fault deformation (Figure 12) shows that there are higher levels in areas with high strain rates, for example, in the SWL and ECSZ.Much of the imaged deformation is in zones that are affected by time-variable processes such as the Long Valley Caldera and Coso magmatic system (Montgomery-Brown et al., 2015;Wicks et al., 2001), and the southern Sierra Nevada, which experiences hydrological loading (Argus et al., 2017;Hammond et al., 2016) and aquifer related deformation (Argus et al., 2017;Hammond et al., 2016;Neely et al., 2021).Non-tectonic processes contribute to the off-fault deformation field because they produce GPS velocity gradients that cannot be well explained as slip deficit on fault systems.
We also see higher values where fault systems terminate and do not continue along another fault system.For example, in the SWL there is a 36 km gap between the southwest end of the Deep Springs Fault and the northern end of the Hunter Mountain/Saline Valley fault system.Here a high amount of off-fault deformation is found in our model because the block models bridge gaps with a boundary that accommodates relative motion between the NSHM faults.In east and northwest Nevada some light-yellow patches indicate deformation inside wide gaps between faults (Figure 12).Better characterizing off-fault deformation is a part of the recommendations of Johnson et al. (2024), and recognizing where the slip rate models fail to explain deformation is the beginning of better accounting for it in future models.

Geodetic Slip Rates Versus Strain Rate
Slip rates tend to be higher in areas with higher strain rates.We can confirm that this is the case in our model by plotting the magnitude of the slip rate as a function of the magnitude of the strain rate near the fault.We define magnitude of strain rate as the norm of the tensor components Kreemer et al., 2014).We exclude faults in Long Valley Caldera because the high strain rates are from non-tectonic processes, and faults below 34.8°latitude because they include the San Andreas fault which has a very high slip rate.Strain rates vary by orders of magnitude in the plate boundary, so we plot the values on a log 10 scale for both axes in Figure 13.The correlation between magnitudes of strain rate and slip rate is particularly strong in log 10 units, having r = 0.72 (for linear scale r = 0.50).
If there is a strong correlation between strain rates and slip rates, this suggests that initial estimates of fault slip rates can be made directly from the strain rate maps.This approach could have utility in initializing deformation models that need prior values for slip rate inversions.These initial values would be independent of geologic slip rates, preserving the independence of the geodetic models.However, we observe that there are many cases where low slip rate faults exist in high strain rate areas.In Figure 13 there are many faults where the strain rates are over 50 × 10 9 /yr but slip rates are below 1 mm/yr or even 0.1 mm/yr.These faults are near other faults that are doing more to accommodate the deformation budget, possibly they are more favorably located, oriented, or are more mature and mechanically efficient.
Conversely, there are no faults in very low strain rate areas that have high slip rates (see the upper left corner of Figure 13).For example, when the strain rate is below 4 × 10 9 /yr there are no faults with slip rates above 0.1 mm/ yr, and when strain rates are below 10 × 10 9 /yr there are no slip rates above 0.6 mm/yr.Thus, the relationship between strain rates and model slip rates suggests that there is an upper bound on slip rate that is a strong function of strain rate.
If we consider the maximum value for slip rate magnitude (s mag ) inside 4 strain rate bins and use them to obtain an line fitting log 10 (ε mag ) versus log 10 (s mag ) we get the red dashed line in Figure 13.This line defines an envelope of the data above which there are almost no values for slip rate magnitude.The formula for the line is: where strain rate magnitude ε mag is in units of 10 9 /yr and slip rate magnitude s mag is in mm/yr.Equation 5 gives, for example, when ε mag = 10 × 10 9 /yr a maximum slip rate bound of 1.12 mm/yr, and when ε mag = 100 × 10 9 the maximum slip rate bound of 6.07 mm/yr.
The red line in Figure 13 describes an envelope based on maximum values of slip rates.However, it is sensitive to outliers since it is based on least squares fit to the maximum slip rates within bins.If we instead take the 95th percentile slip rate inside each of the 4 strain rate bins we get the blue dashed line in Figure 13, whose formula is: s mag,95 < 0.0347 ε mag 0.9922 (6) This gives, for example, when ε mag = 10 × 10 9 /yr then 95% of the slip rates will be below 0.34 mm/yr, and when ε mag = 100 × 10 9 then 95% of the slip rates will be below 3.35 mm/yr.It may be prudent to use the 95th percentile version of the formula since its inference is less sensitive to outlier slip rates.However, it will result in fewer high slip rates which could exist if some faults slip faster than others with similar strain rates for for example, mechanical reasons.Whether there is a similar lower bound on slip rates is equivocal since the lower right area in Figure 13 has faults that slip between 0.001 and 0.1 mm/yr which are very slow and are the most difficult to resolve geodetically.
There are caveats.For example, the slip rates in this study were derived using a regularization that included applying a damping of slip rates that is a function of strain rate, and so could affect the relationships in Equations 5 and 6.However, the bulk of slip rates (∼80%) agree with geological slip rates, suggesting that they are appropriately regularized.The resolution of strain rate maps varies with technique and are sensitive to assumptions about spatial smoothing (Figure S5 in Supporting Information S1) and could lead to differences in Figure 13.Scatterplot of slip rate magnitude versus strain rate magnitude using log 10 scale for both axes.Strain rates at Long Valley Caldera and all faults south of 34.7°are excluded.Red dashed line is slip rate envelope based on maximum of data, blue dashed line is similar but using the 95th percentile of the slip rates within strain rate bins.See Equations 5 and 6 for formulas.
results when using these formulas.The bounding envelopes are tuned to the WL where fault density is relatively high, but the coefficients could be customized for other regions within plate boundary zones where strain rates are higher, more variable or fault density is lower.However, the simplicity of these relationships suggests that when strain rate maps are available, they may be conveniently used to generate bounding a priori values for models of slip rates on faults in complex fault zones.Alternatively, these relations could be used as a check on models, for example, to identify model slip rates that are unusually large outliers in a way that does not rely on geological slip rates.

Conclusions
We have presented a new method for estimating fault slip rates in areas of active tectonic deformation that have many faults that comprise complex networks that may not completely connect to describe contiguous blocks.The method uses an iterative algorithmically driven construction of model block geometries to obtain many geodetic slip rates estimates for each fault.The result obtains better sampling of the epistemic uncertainty associated with limited knowledge of fault connectivity.The method includes constraining the models with a median filtered and interpolated version of the GPS velocity field, applying regularization based on the background strain rate estimated from GPS data.
We applied the method to the WL in the western GB to obtain a robust set of slip rates that agree with geologic slip rates in the USGS NSHM database to within uncertainties ∼80% of the time.This is achieved without constraining the slip rates to be the same as, or in the range of the geologic rates.We also estimated the distribution of off-fault deformation which tends to occur in areas with higher strain rates, areas where faults in the database do not connect end to end to other faults, and in places where non-tectonic signals are present in the GPS velocity data.
Capacity for automatic block model generation and estimating slip rates in a robust way paves the way for larger scale application (e.g., entire Western US or world).Improving slip rate estimates based on geodesy alone will help objectify and strengthen seismic hazards estimates in complex fault systems.

Figure 1 .
Figure 1.(a) Walker Lane and western Great Basin Region of the western United States, red box indicates area covered in panels (b, c), red horizontal dashed lines a-f indicate location of profiles shown in Figure 2. (b) Median spatial filtered velocity field.Vectors are of constant length with color indicating magnitude of velocity in a North America reference frame.(c) Fault networks representing sources in the USGS National Seismic Hazard Maps (Hatem et al., 2022a) (red) and other Quaternary faults (black).Figure S1 in Supporting Information S1 is a version of C annotated with fault names referred to in the text.SNGV indicates the Sierra Nevada/Great Valley microplate.
Figure 1.(a) Walker Lane and western Great Basin Region of the western United States, red box indicates area covered in panels (b, c), red horizontal dashed lines a-f indicate location of profiles shown in Figure 2. (b) Median spatial filtered velocity field.Vectors are of constant length with color indicating magnitude of velocity in a North America reference frame.(c) Fault networks representing sources in the USGS National Seismic Hazard Maps (Hatem et al., 2022a) (red) and other Quaternary faults (black).Figure S1 in Supporting Information S1 is a version of C annotated with fault names referred to in the text.SNGV indicates the Sierra Nevada/Great Valley microplate.

Figure 2 .
Figure 2. Velocity components perpendicular to the profile locations shown in Figure 1.Original velocities are shown with open circles and error bars with two times the uncertainties.Larger blue circles are the median spatial filtered velocities, the small gray dots are values from the gridded velocity field.Velocities parallel to the profile direction are provided in the Figure S4 in Supporting Information S1.Some geographic features are annotated.

Figure 3 .
Figure 3. (a) Shear strain rates ( ė 1 ė 2 ), color scale is in nanostrains (10 9 ) per year.White "C N S B" indicates strain rate anomaly associated with Central Nevada Seismic Belt (CNSB) earthquakes.(b) Same map except with correction applied for viscoelastic postseismic relaxation from all the CNSB earthquakes modeled by Hammond et al. (2009).Black lines are Quaternary faults.Approximate extent of Southern (Southern Walker Lane), Central (CWL) and Northern Walker Lane (NWL) according to Faulds and Henry (2008) are indicated with white curly braces.Fault names are given in Figure S1 in Supporting Information S1.

Figure 5 .
Figure 5. (a) Faults of the southern Walker Lane in the National Seismic Hazard Maps database (blue), other Quaternary faults and CA/NV state line in gray.(b) Initial Delaunay triangulation of nodes that represents primitive block model, (c) blocks after reduction of number of blocks through iterative combination of neighboring blocks, (d) solution for block motion from this model.Color indicates vertical axis spin rate component of the solution Euler pole for each block.Block movement is massively exaggerated to illustrate sense of relative motion, and strain accumulation at block boundaries is removed to emphasis rigid long-term component of motion.

Figure 6 .
Figure 6.(a) Strike slip component of fault slip rates in the Walker Lane (WL), (b) dip slip component (projected to horizontal).Gray colored faults are those for which the National Seismic Hazard Maps database prescribed it to be a vertical strike slip fault, so no dip slip component was estimated.SNGV is the Sierra Nevada/Great Valley microplate.NWL, CWL, and Southern Walker Lane (SWL) denote Northern, Central and SWL.A version of (b) with color scale limits narrowed to 0.1 to 0.1 mm/yr shows variability of the slow slip rates in the Basin and Range of central to eastern Nevada, is presented in the Supplementary Materials.

Figure 7 .
Figure 7. (a) Owens Valley (red) and Independence (blue) Faults.Green circles are locations of nearby GPS stations, gray lines are Quaternary faults.(b) Histogram of estimates of strike slip rate for Owens Valley Fault.No normal rate was estimated since it is categorized as a vertical strike slip fault in the geologic database.The green/ red dashed vertical lines indicate the mean/median values in the distribution respectively.The Independence Fault strike slip rates are shown in panel (c) normal component rates are shown in panel (d).In B, C and D the green/red dotted vertical lines are ±1 standard deviation from the mean and ±1.4826 times the median absolute deviation from the median respectively.

Figure 8 .
Figure 8. Image of median block vertical axis rotation rate.Color scale indicates rate and sign, positive (red) is counterclockwise, negative (blue) is clockwise.

Figure 9 .
Figure 9. Histograms of residual velocity in the east coordinate (left) and north coordinate (right).Top row indicates residuals in mm/yr, bottom row shows residual velocity components normalized by their individual component velocity uncertainties.See text for discussion.

Figure 10 .
Figure 10.(a) Systematicity of residual values S i .Color scale shows blue for randomly oriented residuals and red is for neighboring stations with very similar residual azimuth.(b) is systematicity S' i which normalizes the residual vectors by velocity uncertainty rather than residual vector magnitude.Both measures are unitless.When normalizing using velocity uncertainties the zones of systematic misfit are less extensive and occur areas with significant time-variable deformation.

Figure 11 .
Figure 11.Comparison of geodetic slip rates from this study to (a) slip rates from the geologic slip rate database (Hatem et al., 2022b) excluding the Garlock Fault, (b) same as panel (a) except with log 10 scale axes.The horizontal error bar gives the range of low to high geologic rate.If there is a preferred geologic rate it is plotted with a gray-filled circle, else an open circle is plotted at the mean between the low and high rates.Middle row shows comparison of geodetic slip rates from this study to (c) geodetic slip rates from the National Seismic Hazard Maps (NSHM) geodetic deformation modelers as tabulated by Johnson et al. (2024).(d) is same but on log 10 scaled axes.On all plots the vertical error bars are two times the uncertainty in our geodetic slip rates.(e) shows the geodetic slip rates from this study (green circles), NSHM geodetic rates (red circles), and NSHM geologic rates (gray circles) sorted by slip rate.These curves show the similarity in distribution of the NSHM geologic and NSHM geodetic rates, and more continuous variation of slip rates from this study.Faults with the greatest disagreement between geologic, our geodetic, and NSHM geodetic are annotated in panels (a, c).

Figure 12 .
Figure 12.Map of off-fault deformation rate.Dark green lines are faults from the National Seismic Hazard Maps database (Hatem et al., 2022a), other thin black lines are other faults in the USGS QFFD (USGS and AGS, 2011; USGS and CGS, 2011; USGS and NBMG, 2011; USGS and UGS, 2011).