A New Method to Invert for Interseismic Deep Slip Along Closely Spaced Faults Using Surface Velocities and Subsurface Stressing‐Rate Tensors

Inversions of interseismic geodetic surface velocities often cannot uniquely resolve the three‐dimensional slip‐rate distribution along closely spaced faults. Microseismic focal mechanisms reveal stress information at depth and may provide additional constraints for inversions that estimate slip rates. Here, we present a new inverse approach that utilizes both surface velocities and subsurface stressing‐rate tensors to constrain interseismic slip rates and activity of closely spaced faults. We assess the ability of the inverse approach to recover slip rate distributions from stressing‐rate tensors and surface velocities generated by two forward models: (a) a single strike‐slip fault model and (b) a complex southern San Andreas fault system (SAFS) model. The single fault model inversions reveal that a sparse array of regularly spaced stressing‐rate tensors can recover the forward model slip distribution better than surface velocity inversions alone. Because focal mechanism inversions currently provide normalized deviatoric stress tensors, we perform inversions for slip rate using full, deviatoric or normalized deviatoric forward‐model‐generated stressing‐rate tensors to assess the impact of removing stress magnitude from the constraining data. All the inversions, except for those that use normalized deviatoric stressing‐rate tensors, recover the forward model slip‐rate distribution well, even for the SAFS model. Jointly inverting stressing rate and velocity data best recovers the forward model slip‐rate distribution and may improve estimates of interseismic deep slip rates in regions of complex faulting, such as the southern SAFS; however, successful inversions of crustal data will require methods to estimate stressing‐rate magnitudes.


10.1029/2023EA003069
2 of 15 of a single, planar strike-slip fault (Figure 3a), we determine the spacing of stressing-rate tensors that minimizes the inverse model misfit to the forward model applied slip rate distribution.To assess how well individual and joint inversions of surface velocities and subsurface stressing-rate tensors recover slip along closely spaced and branching faults, we utilize a complex, geologically constrained fault model that simulates the southern SAFS and San Jacinto fault system (SJFS) through the San Gorgonio Pass region (Figure 3b).The SAFS consists of two subparallel pathways for earthquake rupture through the San Gorgonio Pass region, but the relative activity of the two pathways remains a topic of debate (e.g., Blisniuk et al., 2021;Kendrick et al., 2015).Because these two pathways are less than one locking depth apart from one another, inversions of global navigation satellite system (GNSS) velocities alone may not uniquely recover slip-rate distributions along the pathways and at the fault branches.For the complex fault inversion, we intentionally include fault surfaces that are inactive in the forward models to assess how well the inversions can recover zero slip along inactive fault surfaces.The method we present here requires estimates of crustal stressing-rate tensors, which are not currently available.However, the inverse method provides a new approach that may constrain the relative activity of closely spaced parallel faults, such as the two pathways for earthquake rupture through the San Gorgonio Pass, in future applications when crustal stressing-rate tensors are available.

Crustal Data Processing
We utilize focal mechanism-derived stress states and GNSS estimated velocities in southern California for multiple purposes.Previous studies show that long-term forward mechanical models of the SAFS produce slip rates that fit geologic slip rate estimates well (e.g., Cooke & Dair, 2011;Devine et al., 2022;Hatch et al., 2023), and that the interseismic forward model-generated surface velocities agree well with GNSS velocities (e.g., Herbert et al., 2014).Previous studies have not compared stress states generated by a complex SAFS model to focal mechanism-derived stress states.Here, we compare the horizontal maximum compression orientations from interseismic forward models to focal mechanism-derived orientations to further validate a complex SAFS model.Additionally, we use the locations of microseismicity and GNSS stations to assess how deviations from the optimal spacing of data impact the inversions.We also use the data uncertainties to weight the constraining data within the inversions.As the purpose of this study is to test the new approach and stressing-rate tensors are not currently available from crustal data, we do not directly invert the actual GNSS estimated velocities or the focal mechanism-derived stress data, but instead use model-generated data.

GNSS Surface Velocity Locations
We generate surface velocities within the complex SAFS forward models at the locations of 201 permanent GNSS station locations (Figure 1) in the Southern California Earthquake Center's Community Geodetic Model version 1 (Sandwell et al., 2016).We only use the horizontal velocities to constrain the inverse models because this is what would be typically used in GNSS inversions (e.g., Zeng, 2022).

Focal Mechanism-Derived Stress States
Prior to deriving stress information from focal mechanisms of microseismicity, we assess the completeness of and decluster the focal mechanism catalog to remove aftershocks and isolate background seismicity (details provided in Figure S1 in Supporting Information S1; Abolfathian et al., 2019;Martínez-Garzón et al., 2016).We start with 41,110 focal mechanisms from the Southern California Earthquake Data Center from 1981 to 2020 (Hauksson et al., 2012;Yang et al., 2012) that have a nodal plane uncertainty of <45°.Removing focal mechanisms with magnitudes below the limit of completeness reduces bias of small events that occur close to seismic stations but are not represented across the entire region of interest.Following Cooke and Beyer (2018), we calculate the completeness magnitude using the maximum curvature method (Wiemer & Wyss, 2000) and identify three periods with completeness magnitudes that decrease as the density of seismic stations increases.For 1981-2001 the completeness magnitude is 2.0, which decreases to 1.5 for 2002-2011 and to 1.1 for 2012-2020.
To decluster the focal mechanism catalog, we follow the nearest-neighbor approach described by Zaliapin andBen-Zion (2013a, 2013b) and define a nearest-neighbor distance threshold in the space-time-magnitude domain by assessing the distribution of the nearest-neighbor distance for all the events.We exclude events that have a nearest-neighbor distance smaller than the threshold because they may reflect short-term perturbations in the stress field resulting from large events rather than background seismicity.The declustered catalog consists of 10,758 events that have an average fault plane uncertainty of 27° ± 9°.The consistent average slip sense over the 40-year catalog and the consistent rate of seismicity over each completeness magnitude period (Figure S1 in Supporting Information S1) confirms that the declustered catalog represents background seismicity and does not include temporal stress state variations.
The MSATSI code, which is based on the SATSI algorithm (Hardebeck & Michael, 2006), performs formal stress inversions to derive normalized deviatoric stress tensors from groups of focal mechanisms (Martínez-Garzón   -Garzón et al., 2016).The 40-year catalog along the southern SAFS and SJFS yields 54 groups of focal mechanisms from which we derive stress states.Each group consists of focal mechanisms within 7.5 km of a specific location.From 1,000 bootstrap resamplings of the fault plane, we estimate ±10° uncertainty of the orientation of the principal stress axes and 25% uncertainty of the deviatoric stress tensor components.We compare the horizontal maximum stress orientations for the 54 stress states to those of the forward model.

Forward Models
We utilize the Boundary Element Method code Poly3D (Thomas, 1993), which solves the governing equations of continuum mechanics to calculate displacements and stresses within the model to simulate faulting within the crust (e.g., Crouch & Starfield, 1990).The forward models simulate both long-term and interseismic loading of (a) a simple, isolated and vertical strike-slip fault and (b) the complex southern SAFS and SJFS in the San Gorgonio Pass region within a homogeneous and linear-elastic half space (Figure 3).For the complex fault forward models, we utilize the inactive northern slip pathway geometry from Hatch et al. (2023), which is primarily based on the Southern California Earthquake Center's Community Fault Model version 5.3 (Marshall et al., 2021) with   some modifications that improve the model fit to geologic slip rates and uplift (e.g., Fattaruso et al., 2014;Hatch et al., 2023;Herbert & Cooke, 2012).We discretize the fault surfaces into triangular elements that can capture fault curvature and branching.Within all forward models, we prescribe zero opening/closing along all faults.Faults in the long-term forward models intersect a horizontal basal crack at 35 km depth that simulates distributed deformation below the seismogenic zone (Figures S2 and S3 in Supporting Information S1;Marshall et al., 2009).
We simulate interseismic deformation in a two-step back-slip-like approach following Marshall et al. (2009).In the first step, a suite of forward models simulates deformation over several earthquake cycles.Shear-tractionfree faults slip freely in response to loading along far-field horizontal basal patches and slip along nearby faults.The zero shear traction condition simulates low dynamic strength conditions, which is when most of the fault slip occurs (e.g., Di Toro et al., 2006;Goldsby & Tullis, 2011).Following Beyer et al. (2018), we implement an iterative technique to prescribe the desired loading velocity at the model edges (Figure 3).To prevent fault slip rates from artificially going to zero at the lateral edges of the model, we apply slip to driving patches for all faults that extend past the bounds of both models.For the simple fault model of an idealized strike-slip fault, we prescribe far-field loading along the basal crack and apply slip to driving patches that produces a nearly uniform strike-slip rate of 1 mm/year along the vertical fault (Figure 3a).For the complex fault model, we prescribe slip along far-field basal patches consistent with 42 mm/year of far-field loading at an orientation of 322° following Herbert and Cooke (2012) (Figure 3b).Following Beyer et al. (2018), we apply slip rates to driving patches in the complex fault model based on published slip rate estimates for each fault segment (e.g., Fay & Humphreys, 2005;McPhillips & Scharer, 2018;Meade & Hager, 2005;Sharp, 1981;Weldon & Sieh, 1985).
In the second suite of forward models, we apply the long-term model slip rates below a prescribed locking depth to simulate interseismic deformation.For the simple fault model, we test the inverse approach with forward model locking depths of 10, 15, and 20 km.For the complex fault model, we utilize a locking depth of 20 km based on the maximum depth of seismicity across the San Gorgonio Pass region (e.g., Yule & Sieh, 2003).To reduce artifacts that would result from an abrupt change in prescribed slip at the locking depth, we create a transition zone by prescribing half of the long-term slip rate to elements that have centroids within 2.5 km of the locking depth.This study tests if the new inverse approach can recover deep interseismic slip rates along complex fault geometries that include closely spaced and branched faults.For simplicity of this test, the complex interseismic model only applies deep slip from the first suite of forward models along the primary faults in the region, the San Andreas and San Jacinto faults.The interseismic models produce surface velocities and stressing-rate tensors at regularly spaced points for both the simple and complex fault models.Within the complex fault model, we additionally query surface velocities at specific GNSS station locations and stressing-rate tensors at locations of recorded microseismicity.To compare the interseismic principal stress orientations with those derived from crustal focal mechanisms, the model includes all of the faults shown in Figure 1, not only the primary faults.

Inverse Models
We use the MATLAB code triinv (Loveless, 2023), which is based on algorithms from Meade (2007), to calculate partial derivatives that relate the stressing rates or surface velocities at specific locations to unit slip rate on each triangular dislocation element within each model.Because MSATSI produces normalized deviatoric stress tensors, we set up separate inversions for the forward model-generated full, deviatoric, and normalized deviatoric stressing-rate tensors.For deviatoric and normalized deviatoric stressing-rate tensor inversions, we remove the mean stress component of the partial derivative.Laplacian smoothing within the inversions prevents abrupt steps in slip rates that would not be expected along crustal faults.We test a range of smoothing weighting parameters to optimize the surface velocity, stressing rate, and joint inverse model performance.The results of the smoothing parameter value testing are independent of the surface velocity and stressing-rate tensor spacing.Within all inversions, elements in direct contact with the free surface of the model (0 km depth) are locked and opening/closing is prohibited.However, we do not constrain the locking depth or sense of slip on any faults in the inverse models.
We assess the performance of individual and joint inversions that use forward model-generated surface velocities and stressing-rate tensors.The simple fault model allows us to determine the optimal stressing-rate tensor configuration and smoothing weight.Inversions of regularly gridded surface velocities have 10 km spacing, which is based on the approximate current permanent GNSS station density in the San Gorgonio Pass region (Figure 1).We test 60 stressing-rate tensor configurations that are based on the microseismicity in the San Gorgonio Pass region, which generally occurs above 20 km depth.Because each stressing-rate tensor represents a potential centroid of a group of microseismic focal mechanisms with a radius between 2.5 and 7.5 km, we limit the stressing-rate tensor depths to between 15 and 7.5 km.All the stressing-rate tensor configurations include either a single row of tensors at a single depth (7.5, 10, 12.5, or 15 km) or two rows of tensors at two separate depths (7.5 and 15 km) on either side of the simple fault.To reduce overlap of focal mechanisms within each group, we define a 10 km minimum along-strike spacing of stressing-rate tensors and only test two rows for stressing-rate tensors at 7.5 and 15 km depths.To reduce the chance that a focal mechanism group would include microseismicity on both sides of the same fault, all stressing-rate tensor locations are at least 5 km away from the fault.We assess the same spacings for the simple interseismic forward model with three different locking depths: 10, 15, or 20 km; this allows us to assess the impact of locking depth on the stressing-rate tensor configuration that best recovers the forward model slip rates.
We use the complex fault model to assess the performance of inversions on a geometrically complicated fault system consisting of multiple closely spaced (<12 km) and interconnected faults.We invert the forward model-generated stressing-rate tensors and surface velocities using a model with two slip pathways from Hatch et al. (2023) to assess how well the inversions recover slip along the portion of the northern slip pathway that is inactive in the forward models.The complex fault model inversions utilize regularly spaced surface velocities and the configuration of stressing-rate tensors that optimizes the simple fault model inversion performance as well as surface velocities at GNSS station locations and stressing-rate tensors at locations of microseismicity groups.We prescribe an uncertainty of 0.3 mm/year to all surface velocity components, which is based on the lowest estimates of GNSS errors for stations that we include (Sandwell et al., 2016).We query stressing-rate tensors at 100 locations following the optimal distribution informed by the simple fault model.Inverse models utilize either all 100 tensors or only 54 tensors at locations with a minimum of 40 nearby cataloged focal mechanisms, which allows for a robust stress state estimate.We prescribe a conservative uncertainty of 25% to all stressing-rate tensor components, at the high end of the estimated uncertainty.When describing the inversions that use only the 54 stressing-rate tensors at locations with a minimum of 40 nearby focal mechanisms and the surface velocities at locations of GNSS stations, we refer to these inversions as using crustal limited locations or as crustal limited inversions.
To assess how well each inversion of forward model-generated stress rate and velocity predictions recovers the prescribed fault slip rates, we calculate the misfit of the inverse model slip rate distribution to the forward model applied slip rate.Because the root-mean-square error can overestimate the model error by emphasizing outliers (Willmott et al., 2017), we define the model performance based on the inverse model misfit to the forward model slip distribution with the area-weighted average misfit per element using Equation 1, where j is the number of elements, S I is the inversion estimated slip rate for an element, S F is the forward model slip rate for an element, and A is the area for an element.

Determination of the Optimal Stressing-Rate Tensor Spacing
An assessment of 60 different stressing-rate tensor configurations reveals the spatial configuration of stressing-rate tensors that best recover the forward-model slip rate distribution (Figures 4a-4e).Figures 4 and 5 present results from inversions of stressing-rate tensors and surface velocities generated by a forward model with a 15 km locking depth, and the Supporting Information S1 (Figure S4 in Supporting Information S1) contains results from the models with 10 and 20 km locking depths.The forward-model prescribed locking depth does not significantly impact the optimal stressing-rate tensor spacing (Figure S4 in Supporting Information S1).Twenty-three of the 60 stressing-rate tensor spacings that we test produce misfits less than or equal to the surface velocity inversion misfit of 0.08 mm/year.Increasing the tensor depth and distance from the fault generally improves the inverse model performance (Figures 4a-4d).Inverting stressing-rate tensors at two separate depths rather than at a single depth improves model performance (Figures 4a-4e).Inverting stressing-rate tensors at both 7.5 and 15 km depth at points that are 5 or 10 km away from the fault with along-strike spacing of 10 km best recover the forward model prescribed slip rate distribution (Figures 4a-4e and Figure S4 in Supporting Information S1).As the along-strike spacing increases to 15 and 20 km, the inverse model performance generally decreases.
We present the smoothing parameter value assessment results from inversions that utilize two rows of stressing-rate tensors at 7.5 and 15 km depths that are 10 km away from the fault with 10 km along-strike spacing (Figure 4f).
Varying the smoothing parameter impacts both the inversion misfit and condition number.A lower condition number indicates the inversion has greater numerical stability.Because using a smoothing parameter value of 0.1 produces misfits within 2% of the minimum misfit and a condition number three orders of magnitude lower than the inversions that produce minimum misfits (Figure 4f), we use this smoothing parameter value for all the inversions.

Assessment of the Inversion Performance
We compare the area-weighted average element misfit for the portion of the fault displayed in Figure 5a to determine which inverse model best recovers the forward model slip rate distribution (Figure 5b).The inversions that use surface velocities and stressing-rate tensors that include magnitude recover both the magnitude and pattern of forward model slip rates well (Figures 5c-5g).Even without prescribing a locking depth within the inversion, the inverse models recover the forward-model locking depth well.The inversions estimate a broader locking depth transition zone than is prescribed in the forward model, but the inversions recover slip rates slower than 0.1 mm/year for all elements above 10 km, which are locked in the forward model (Figure 5).The inversion of the surface velocities produces a misfit of 0.08 mm/year, which exceeds that of the stressing-rate tensor inversion of 0.06 mm/year.The joint inversion that utilizes both full stressing-rate tensors and surface velocities outperforms both individual inversions producing a misfit of 0.04 mm/year.
The largest difference between the inverse models and the forward model applied slip rates are along elements with at least one vertex at the locking depth of 15 km (Figures 5c-5h).The inversions overestimate slip on elements just above the locking depth transition zone and underestimate slip on elements within and below the locking depth transition zone.This result highlights the limit of this inverse approach to capture sharp changes in slip rate along faults due to the applied Laplacian smoothing.Because we do not have evidence that locking depth transition zones within the crust are as sharp as we prescribe in the forward models, this smoothing across the locking depth does not cause concern.However, implementing a sparsity-promoting regularization instead of Laplacian smoothing could better recover sharp changes in slip rates (e.g., Evans & Meade, 2012).
Because current methods of deriving stress information from focal mechanisms produce normalized deviatoric stress tensors (e.g., Martínez-Garzón et al., 2014), we assess the performance of inverse models that use either deviatoric or normalized deviatoric stressing-rate tensors.These inversions reveal the impact of removing the mean normal stress component and stress magnitude from the inverse model constraint.Removing the mean normal stress from the full stressing-rate tensor does not significantly impact the inverse model performance.The deviatoric stressing-rate tensor inversion produces a misfit equal to that of the full stressing-rate tensor inversion (0.06 mm/year).The normalized deviatoric stressing-rate tensors lack magnitude; thus, the inversion is poorly posed to recover slip rates with magnitude.As we expect, removing the stressing-rate tensor magnitude leads to the inverse model estimating near zero slip rates along the entire fault (Figure S5 in Supporting Information S1).
Consequently, the inversion recovers the locked, shallow portion of the fault well but not the deep slip rates or the locking depth.Because the normalized deviatoric stressing-rate tensor inversion for the simple fault model failed to recover the forward model slip rate distribution, henceforth, we only discuss results from model inversions that use full or deviatoric stressing-rate tensors that include magnitude.
Overall, the joint inversions recover the forward model slip better than or as well as the individual inversions (Figure 5).Although the individual deviatoric and full stressing-rate tensor inversions perform similarly, the joint inversion that utilizes the deviatoric stressing-rate tensors does not recover the slip rates near the locking depth transition zone as well as the joint inversion that utilizes the full stressing-rate tensors.Still, both joint inversions recover the forward model slip rate distribution better than or as well as all the individual inversions.

Forward Model Validation
To validate the complex forward fault models, we compare the maximum horizontal compression orientation for the model and focal mechanism-derived stress tensors (Figure 6).At 29 of the 54 crustal locations, the forward interseismic model produces maximum horizontal compression orientations that are within 2 standard deviations (3°-15°) of the crustal orientations.The stress states derived from focal mechanisms show spatial variations in the maximum horizontal compression orientation whereas the forward model-generated stressing-rate tensors produce relatively uniform approximately north-south oriented maximum horizontal compression orientations across the region of interest.Most of the locations where the model results do not match the crustal data well are at 7.5 km depth and near the inactive portion of the northern slip pathway (Figure 6).Where the model results differ from crustal data, the model may not completely capture the crustal faulting behavior.For example, some fault structures may be oversimplified or missing from the model, such as the Cox Ranch and Beaumont Plain fault zones (e.g., Yule & Sieh, 2003), which could impact the maximum horizontal compression orientation at specific locations.Further exploration of the activity and geometry of faults along and near the northern slip pathway along the SAFS in the San Gorgonio pass region may provide insight on how to improve the model fit to the crustal data.Another potential explanation for the differences between the crustal and modeled maximum horizontal compression orientation is that the focal mechanism derived stress states represent the total stress state unlike the modeled stressing-rates that solely capture interseismic stresses.Overall, the forward model results are consistent with regional studies that invert focal mechanisms for the entire area and show approximately north-south oriented horizontal maximum compression (e.g., Hardebeck & Hauksson, 2001).

Inverse Model Results
We present results from inversions of forward model-generated deviatoric stressing-rate tensors and surface velocities that are either regularly spaced or only at locations where data is currently available from the southern California focal mechanism catalog and GNSS stations (Figures 1 and 7).Similar to the simple fault model inversions, all the complex inverse models recover the approximate locking depth applied in the forward model.For all the complex fault model inversions, the area-weighted average element misfit increases with depth until ∼22.5 km depth, below which the average misfits remain high (Figure 8a).In general, the misfit for the joint inversions increases less with depth compared to the individual inversions, meaning that for the joint inversions, the resolution of slip rates is more equal at all depths compared to individual inversions (Figure 8a).As a consequence of the smoothing, the inversion underestimates slip rates below the locking depth.Because this misfit is pervasive across the entire model and is not localized to one fault strand or segment, the overall misfit with depth is generally largest within 5 km of the 20 km locking depth (Figure 8a).The joint inversions produce smaller misfits than both individual inversions that use regularly spaced and crustal limited locations (Figure 8).
To determine which inversion of regularly spaced data best recovers the forward model slip distribution for the entire region of interest, we compare the area-weighted average element slip rate misfit (Equation 1; Figure 8).The regularly spaced surface velocity inversion produces an overall slip rate misfit of 1.4 mm/year, which is slightly larger than the 1.3 mm/year misfit of the regularly spaced deviatoric stressing-rate tensor inversion (Figure 8b).The regularly spaced stressing-rate tensor inversion recovers forward model slip better above and within the locking depth transition zone than the regularly spaced surface velocity inversions (Figure 8b).Inverting the regularly spaced data jointly produces the lowest misfit (1.0 mm/year; Figure 8b).
Inversions that utilize stressing-rate and velocity data only at crustal limited locations generally recover the forward model locking depth and slip rate distribution (Figure 8).For individual inversions, inverting crustal limited deviatoric stressing-rate tensors produces a larger misfit than the crustal limited surface velocity misfit (1.8 > 1.4 mm/year).Below the locking depth, the inversion of deviatoric stressing-rate tensors at crustal limited locations does not recover deep slip rates as well as the inversion of surface velocities at GNSS station locations (Figure 8).The crustal limited joint inversion produces a lower misfit (1.2 mm/year) than the individual crustal limited and regularly spaced inversions (Figure 8b).
Inverting regularly spaced stressing-rate tensors and surface velocities improves the overall inversion performance compared to inverting only information at crustal limited locations.The regularly spaced surface velocity inversion includes 198 surface velocity locations, and the crustal limited surface velocity inversion includes 10.1029/2023EA003069 10 of 15 201 locations.The small difference in the number of constraining data may explain the similar misfit of both surface velocity inversions, but the difference in spatial distribution of the constraining data could contribute to the differences in the misfits along individual fault strands or segments (Figure 8).Reducing the number of deviatoric stressing-rate tensors that constrain the individual inversions from 100 to 54 leads to an overall increase in the inverse model misfit to the forward model slip-rate distribution.Furthermore, the 54 deviatoric stressing-rate tensor crustal limited locations are not evenly distributed across the region of interest.A significant gap in microseismicity along the southern SAFS reduces the number of stressing-rate tensors constraining the inversion by 33% (Figures 1 and 6).This reduction could explain why the crustal limited deviatoric stressing-rate tensor inversion cannot resolve slip rates along some fault segments as well as the regularly spaced deviatoric stressing-rate tensor inversion.
We expect the largest misfits around fault branches and along closely spaced faults where inversions cannot uniquely resolve slip rates.The San Bernardino segment directly connects to both the inactive portion of the northern slip pathway and the active southern pathway of the southern SAFS forming a branched fault (Figure 1).Comparing the slip rate misfits along individual fault segments and strands provides insight on how well each inversion can recover slip rates at fault branches and along the two subparallel slip pathways of the southern SAFS.The San Bernardino segment of the SAFS yields the greatest misfit for all the inversions (Figure 8b).Due to smoothing of slip rate across faults within the inversion, the inverse models overestimate slip rates along the inactive portion of the northern pathway (Figure 9 red colors) and underestimate slip rates along the adjacent San Bernardino segment (Figure 9 blue colors).The tradeoff in slip rates among the branched fault segments is lesser for the joint inversion.As a result, the joint inversion misfits along the inactive portion of the northern pathway and the San Bernardino segment are smaller than the misfits for the inversions of individual constraints.

Constraint Weighting in Joint Inversions
The weighting of surface velocities and stressing-rate tensors within the inversions depends on three parameters: (a) the relative numbers of constraint components, (b) the prescribed uncertainties, and (c) smoothing weighting.
Because multiple factors impact the weighting of differing data types, the surface velocities and stressing-rate tensors are likely not equally weighted in the joint inversions.Each stressing-rate tensor consists of six components (three shear and three normal), and each surface velocity consists of two components (east and north).
For the regularly spaced joint inversions, a greater number of stressing-rate tensor components constrain the  inversion than surface velocity components; this means that the stressing-rate tensors may have more weight in the joint inversion than the surface velocities.In contrast, for the crustal limited joint inversions, a greater number of surface velocity components constrain the inversion than stressing-rate tensor components.Regardless of the ratio of stressing-rate tensor to surface velocity components constraining the inversions, increasing the amount of constraining information improves the inverse model's recovery of forward model slip rates.Increasing the number of surface velocity locations by utilizing campaign GNSS stations or InSAR data could potentially improve the inversion performance.The second factor that impacts the weighting of the two data types is the uncertainty we prescribe to each component.Because each component for surface velocities and stressing-rate tensors has uncertainty of 20%-40% of the component, the two data types have similar weighting in the joint inversions.Since the smoothing weighting can also impact how the inverse model constraining information is weighted, we assess the impact of varying the smoothing weighting on the slip rate misfit for the complex fault inversions.We find that a range of smoothing weightings (varying by a factor of 10 4 ) for all the inversions produce slip rate misfits that vary by <0.05 mm/year (Figure S6 in Supporting Information S1), which suggests that the inversions are more sensitive to the number and location of constraining data than the smoothing weighting.

Comparison of Individual Inverse Model Results
The regularly spaced stressing-rate tensor inversions may have better overall performance than the surface velocity inversions because the stressing-rate tensors are at depth, closer to the locking depth transition zone and the slipping portion of faults.The stressing-rate tensor spacing assessment shows that for inversions that utilize stressing-rate tensors at a single depth the misfit generally decreases as the stressing-rate tensor depth increases.Many of the simple fault model stressing-rate tensor inversions that utilized tensors at a single depth outperformed the surface velocity inversion, and the addition of stressing-rate tensors at a second depth further improved the stressing-rate tensor inversion performance.Furthermore, the joint inversions include constraints at three separate depths (0, 7.5 and 15 km) and best recover forward model slip rates for both the simple fault and complex fault models.More information on spatial variations of conditions may yield more accurate inversions for slip rate.
For the complex fault model, the surface velocity inversions can recover deep interseismic slip rates (>25 km depth) better than stressing-rate tensor inversions (Figure 8a).The assessment of the optimal spacing of stressing-rate tensors shows that decreasing the along-strike tensor spacing from 20 to 10 km can improve the inversion performance (Figures 4a-4e), suggesting that stressing-rate tensors may provide higher resolution slip rate information over short distances (10-15 km).Consequently, the stressing-rate tensors provide better slip rate information along portions of faults closest to the tensors (<25 km depth) than below the locking depth.Even though the interseismic surface velocities are farther from the slipping portions of faults than the subsurface stressing-rate tensors, the ability of the surface velocities to resolve slip rates is less sensitive to their distance from the fault.As a result, surface velocity inversions may better constrain interseismic slip rates along deep portions of the fault (>25 km depth) than stressing-rate tensor inversions (Figure 8a).In addition to having a greater number of inputs, the joint inversion takes advantage of the benefits of both data types, which improves the inverse model performance compared to individual inversions (Figure 8).

Future Application to Natural Fault Systems
The complex fault models show that joint inversions of stressing-rate tensors and surface velocities could improve current estimates of slip rates along closely spaced and branching faults; the distribution of these rates can help constrain both the locking depth and relative activity of closely spaced faults.For example, joint inversions resolve slip rates well along the northern pathway of the southern SAFS through the San Gorgonio Pass where fault activity remains debated (e.g., Blisniuk et al., 2021;Kendrick et al., 2015).
Implementing the inverse method that we present here for any crustal fault system requires a priori information including geodetic and microseismic catalogs as well as a three-dimensional fault geometry, and uncertainty or inaccuracy in the inverse model inputs propagates through the model.Because we invert forward model generated stressing-rate tensors and surface velocities, we know that the fault geometry used in the inversions is accurate.As a consequence, the inversion misfits that we calculate exclude uncertainty that may stem from uncertainty or inaccuracy in the model fault geometry.In addition to uncertainty related to the a priori information, model parameters, such as fault element size, may impact the inverse model performance.In this study, the simple and complex fault models have average element lengths of 3-5 km.Future applications of the inverse method we present here should consider that the average element length could impact the optimal stressing-rate tensor spacing.
Because microseismicity in the crust is generally not evenly distributed across a region (Figure 1), the optimal regular spacing that we determine from the idealized simple fault model may not be available for crustal data sets.For the complex SAFS model, limiting the stressing-rate tensor locations to points with sufficient nearby recorded focal mechanisms increases the average misfit of the joint inversion, but the inversion still estimates <2.0 mm/year of strike-slip along the inactive northern pathway.With time and additional microseismicity, focal mechanism catalogs may enable additional tensor locations to be included in the model, which would improve the spatial consistency in model performance.
Another challenge prevents us from applying this new method to crustal data at this time: we do not know of a method to reliably estimate deviatoric stress magnitude and stressing rate within the crust.The results of this study show that inversions of deviatoric stressing-rate tensors perform similarly to inversions that utilize full stressing-rate tensors, meaning that inversions of crustal data would not require mean normal stress state information.The stress states inferred from focal mechanisms provide normalized stress due to microseismicity but not magnitudes.A recent study provides a method to estimate absolute stress magnitude from focal mechanisms and precisely located earthquakes (Fialko, 2021).However, absolute stress does not directly correspond to interseismic stressing rates that are necessary to invert for slip rates.Absolute stress evolves with time since the last earthquake so that microseismicity responds to the total stress state, which includes the effect of accumulated tectonic loading, not solely stressing rates from interseismic loading.If we can derive crustal deviatoric stressing rates, then we may be able to provide additional constraint on deep slip rates along faults in the San Gorgonio Pass region, which would reveal locking depths and relative fault activity.

Conclusions
We present a new method that utilizes model-generated interseismic surface velocities and subsurface stressing-rate tensors to estimate three-dimensional slip rate distributions along a simple, isolated strike-slip fault model and a complex fault model that simulates the southern SAFS.The inversions of forward model-generated stressing-rate tensors and surface velocities for the simple fault model reveal that a sparse, regularly spaced distribution of stressing-rate tensors can recover the forward model slip rate distribution better than surface velocity inversions alone.Additionally, inversions that utilize deviatoric stressing-rate tensors recover the slip rates along faults as well as inversions of full stressing-rate tensors.Inverting forward-model-generated surface velocities and subsurface stressing-rate tensors jointly recovers both the simple and complex forward model applied slip rate distributions better than inverting velocity and stress information individually.For the complex fault model that simulates the SAFS through the San Gorgonio Pass region, inversions of regularly spaced velocity and stress information recover the forward model slip rates better than inversions of velocity and stress information only at locations where crustal data is currently available.
Future joint inversions of surface velocities from GNSS stations and subsurface deviatoric stressing rates potentially derived from microseismic focal mechanisms could provide additional constraint on the deep slip distribution and as a result both the interseismic locking depth and relative activity of faults along closely spaced faults.The complex fault inversions generally recover very slow slip rates along the northern pathway of the SAFS that is inactive in the forward model, suggesting that the method we present here could be used to inform the activity of the northern and southern pathways of the SAFS through the San Gorgonio Pass.However, prior to applying this new method to invert crustal data sets, we require a method to reliably estimate the deviatoric stressing rates that include magnitude.With an increase in the number of available microseismic focal mechanisms with time and a method to calculate stressing rates from focal mechanisms, the method we present here could improve constraints on fault slip rate distributions in regions with closely spaced and branching faults.

Figure 1 .
Figure 1.Map of the San Gorgonio Pass region with the modeled fault surface traces for the region of interest.Black, orange (active northern pathway), blue (southern pathway), and light blue (San Bernardino segment of the San Andreas fault system) traces indicate active faults in all complex forward models.Red traces (inactive northern pathway) indicate faults that are inactive in all forward models.Gray traces indicate the secondary faults that are active in the long-term forward models only.Teal open circles show microseismicity from the declustered catalog.White triangles show global navigation satellite system stations that we use.White box shows the area we use to calculate the inverse model misfits.

Figure 2 .
Figure 2. Flow chart showing (a) the methods we use to assess the new inverse method and (b) the steps for a future application of the inverse method.Polygons on the left of a model are inputs.Polygons on the right of a model are outputs.Parallelograms indicate a model output is used as an input in the next model.

Figure 3 .
Figure 3.The long-term forward model geometries of the (a) simple and (b) complex fault models.(a) The green surface indicates the area we use to calculate the inverse model misfit and show in Figures 4 and 5b-5h.Rates adjacent to extended fault patches indicate the applied slip rates.Arrows on the far-field basal crack show applied loading.(b) Modified from Beyer et al. (2018).Arrows indicate the applied tectonic velocities along the far-field basal crack.
Figure 3.The long-term forward model geometries of the (a) simple and (b) complex fault models.(a) The green surface indicates the area we use to calculate the inverse model misfit and show in Figures 4 and 5b-5h.Rates adjacent to extended fault patches indicate the applied slip rates.Arrows on the far-field basal crack show applied loading.(b) Modified from Beyer et al. (2018).Arrows indicate the applied tectonic velocities along the far-field basal crack.

Figure 4 .
Figure 4. (a-e) Each square represents one stressing-rate tensor spacing with the color indicating the average element misfit.We invert one row of stressing-rate tensors at (a) 7.5, (b) 10, (c) 12.5, or (d) 15 km depth or two rows of stressing-rate tensors at (e) 7.5 and 15 km depths.The red box indicates the optimal spacing.(f) Average element misfit (left y-axis) and joint inversion condition number (right y-axis) against smoothing parameter for inversions that use surface velocities with 10 km spacing and stressing-rate tensors with the optimal spacing (red box in (e)).The black line shows the minimum misfit for the joint inversion.The gray rectangle indicates the smoothing parameter value we use.

Figure 5 .
Figure 5. (a) The 3-D fault model geometry with the optimal stressing-rate tensor spacing and the grid of surface velocities.(b-h) Show the strike-slip rate or strike-slip rate difference for the patch shown in (a).(b) The 15 km locking depth interseismic forward model applied strike-slip rates.(c-h) The difference between the forward model applied strike-slip rates and the inversion estimated strike-slip rates.Blue indicates the inversion underestimates slip rates and red indicates the inversion overestimates slip rates.

Figure 6 .
Figure 6.Maximum horizontal compression orientation for the focal mechanism derived normalized deviatoric stress tensors (S crust , gray lines) and the forward model generated stressing-rate tensors (S model , green/orange lines) at 7.5 (a) and 15 km (b) depths.Green lines indicate the model results are within 2 standard deviations (std) of the focal mechanism derived results.Circle color shows 2 std of the focal mechanism derived results.Gray dots indicate locations of optimally spaced stressingrate tensors that did not have enough nearby focal mechanisms to estimate the crustal stress tensor.Black fault traces show active faults and red traces show inactive faults in the forward models.

Figure 7 .
Figure 7. Locations of regularly spaced (a) surface velocities (triangles) and (b) stressing-rate tensors (circles) for the complex fault model.Red fault trace indicates the inactive portion of the northern slip pathway in forward models.(a) Map view with black box indicating the region used for misfit calculations.(b) Oblique view of San Andreas fault system and San Jacinto fault system geometry colored by the forward model slip rates.

Figure 8 .
Figure 8.The area-weighted average element misfit (a) with depth and (b) for the entire region of interest and individual fault segments (defined in Figure 1) for the deviatoric stressing-rate tensor (light blue), surface velocity (orange) and joint (indigo) inversions.(a) Each point is the misfit for elements within 2.5 km of the specified depth.Solid lines-regularly spaced inversions.Dashed lines-crustal limited inversions.(b) Vertical lines-regularly spaced inversions.Open circlescrustal limited inversions.

Figure 9 .
Figure 9.The difference between the forward model applied and the regularly spaced stressing-rate tensor inversion estimated strike-slip rates along a portion of the inactive northern pathway (INP; red outline) and the San Bernardino segment (black outline).(a) Map of region of interest with gray box indicating the area shown in perspective views in (b) (from the south) and (c) (from the north).(b, c) Red elements indicate the inverse model overestimates slip rates while blue elements indicate the inverse model underestimates slip rates.Fault elements are transparent above the 20 km locking depth.