Abstract

Present-day shear tractions along faults of the San Gorgonio Pass region (southern California, USA) can be estimated from stressing rates provided by three-dimensional forward crustal deformation models. Due to fault interaction within the model, dextral shear stressing rates on the San Andreas and San Jacinto faults differ from rates resolved from the regional loading. In particular, fault patches with similar orientations and depths on the two faults show different stressing rates. We estimate the present-day, evolved fault tractions along faults of the San Gorgonio Pass region using the time since last earthquake, fault stressing rates (which account for fault interaction), and coseismic models of the impact of recent nearby earthquakes. The evolved tractions differ significantly from the resolved regional tractions, with the largest dextral traction located within the restraining bend comprising the pass, which has not had recent earthquakes, rather than outside of the bend, which is more preferentially oriented under tectonic loading. Evolved fault tractions can provide more accurate initial conditions for dynamic rupture models within regions of complex fault geometry, such as the San Gorgonio Pass region. An analysis of the time needed to accumulate shear tractions that exceed typical earthquake stress drops shows that present-day tractions already exceed 3 MPa along portions of the Banning, Garnet Hill, and Mission Creek strands of the San Andreas fault. This result highlights areas that may be near failure if accumulated tractions equivalent to typical earthquake stress drops precipitate failure.

INTRODUCTION

The southern San Andreas fault system (southern California, USA) consists of multiple active faults that accommodate the deformation between the North American and Pacific plates. Accurate estimates of the earthquake hazard in California require an accurate assessment of the potential for large throughgoing earthquakes and the ability of ruptures to propagate through fault intersections and complexities (e.g., Field et al., 2014). One region of such complexity is the San Gorgonio Pass region (SGPr), the site of a restraining stepover along the southern San Andreas fault (Fig. 1). Accurate dynamic rupture models of the SGPr that simulate potential rupture paths will help us assess the potential for large and damaging earthquakes through this region (e.g., Tarnowski, 2017; Douilly et al., 2020).

Dynamic rupture models show that, in general, the size and extent of earthquake ruptures can depend highly on the initial conditions of the model (e.g., Oglesby, 2005). These conditions include physical aspects, such as fault geometry and location of rupture nucleation (e.g., Lozos et al., 2012; Lozos, 2016; Tarnowski, 2017), and time-dependent aspects, such as state of stress and frictional parameters (e.g., Kame et al., 2003; Aochi and Olsen, 2004; Kase and Day, 2006; Duan and Oglesby, 2007). Dynamic rupture models typically prescribe initial shear and normal tractions by resolving the remote stress tensor, constrained from focal mechanism inversions, onto individual fault elements (e.g., Kame et al., 2003; Oglesby et al., 2003). This approach provides spatially variable “resolved” tractions that capture the first-order loading of the faults, but does not take into account the loading history nor the prior stress interactions between faults. Not only can individual earthquake events change tractions along nearby faults, advancing or retarding each fault’s earthquake clock (e.g., King et al., 1994; Stein, 1999; Duan and Oglesby, 2005), but interaction among neighboring active faults influences their long-term slip rates and stressing rates (Willemse and Pollard, 1998; Maerten et al., 1999; Loveless and Meade, 2011). Interseismic stressing rates on faults can be estimated using geodesy (e.g., Smith and Sandwell, 2006). Smith and Sandwell (2006) calculated the Coulomb stress along the southern San Andreas fault just prior to the A.D. 1812 and 1857 earthquakes, showing elevated Coulomb stress in the predicted epicenter locations for these earthquakes. However, the total accumulated traction along any given fault segment depends on both the accumulated tractions during the interseismic period as well as nearby rupture history (e.g., Smith-Konter and Sandwell, 2009; Richards-Dinger and Dieterich, 2012; Tong et al., 2014). For example, Smith-Konter and Sandwell (2009) found that California fault segments with higher rates of interseismic loading have shorter earthquake recurrence intervals.

To account for both loading history and fault interaction and to produce reasonable estimates of fault stress, we simulate deformation within the San Gorgonio Pass region using three-dimensional forward models that provide both steady-state slip rates over multiple earthquake cycles and stressing rates between earthquake events. Because we use steady-state slip rates over multiple earthquake cycles to drive models that simulate interseismic deformation, the resulting shear stressing rates incorporate the interactions between faults of the southern San Andreas fault system. The interseismic shear stressing rates along with information about time since last earthquake event can be used to estimate the shear traction on faults through the SGPr following the approach employed by Tong et al. (2014). The resulting estimates of shear traction may differ from those resulting from resolving the remote stress tensor onto faults, in that our models explicitly include fault interaction and fault loading from depth during the interseismic period. Furthermore, we incorporate the effects of recent earthquakes on faults near the SGPr to produce a more accurate estimate of the current stress state of this system. Using tractions that incorporate fault interaction and loading history may enhance the accuracy of dynamic rupture models, refining our insight into the nature of potential earthquake rupture propagation within the SGPr.

REGIONAL GEOLOGY

The San Gorgonio Pass Region

Through the San Gorgonio Pass region (SGPr), deformation is partitioned onto multiple active and nonvertical fault strands (e.g., Matti et al., 1992; Fig. 1). The San Bernardino strand of the San Andreas fault lies at the northwestern end of the San Gorgonio Pass. Two potential rupture pathways go through the restraining bend connecting the San Bernardino strand to the Coachella segment of the San Andreas fault. The southern pathway consists of the San Gorgonio Pass thrust, Garnet Hill strand, and Banning strand of the San Andreas fault (Fig. 1). The San Gorgonio Pass thrust dips to the north and has a corrugated geometry near the Earth’s surface (e.g., Matti et al., 1992). The eastern end of the San Gorgonio Pass thrust connects to the Garnet Hill and Banning strands. The north-dipping Garnet Hill and subparallel Banning strands have approximately the same strike. The northern pathway through the SGPr consists of the Mill Creek, Mission Creek, and Galena Peak strands of the San Andreas fault (Fig. 1). Ongoing debate centers on the geometry and activity of these fault strands through the northern part of the SGPr (e.g., Kendrick et al., 2015; Beyer et al., 2018; Fosdick and Blisniuk, 2018). The San Jacinto fault is subparallel to the San Andreas fault and extends to within 2 km of the San Andreas fault at Cajon Pass, northeast of the San Gorgonio Pass. While the San Jacinto fault lies outside of the SGPr, it interacts with the San Andreas fault and consequently impacts both long-term slip rates (e.g., Herbert et al., 2014) and earthquake rupture paths (Lozos, 2016) on the San Andreas fault.

Recent paleoseismic data within the SGPr suggest that previous large throughgoing earthquakes have a recurrence interval of ∼1000 yr, with the most recent earthquake rupture through the San Gorgonio Pass along the southern pathway in ∼1400 A.D. (Heermance and Yule, 2017). Earthquakes along the San Bernardino strand (northwest of the restraining bend) and Coachella segment of the San Andreas fault (southeast of the bend) occur more frequently, with recurrence intervals of 200–300 yr (e.g., Philibosian et al., 2011; Field et al., 2014; Onderdonk et al., 2018). The difference in recurrence intervals outside of and inside of the restraining bend suggests that previous earthquakes that have ruptured along the San Bernardino and Coachella segments terminated at the restraining bend, which may be acting as an “earthquake gate”. During the interseismic period since the last rupture event through the bend, shear tractions have been accumulating along faults within the SGPr. Furthermore, recent earthquakes along faults surrounding the SGPr could have impacted the state of stress within the San Gorgonio Pass, and thus the shear and normal tractions along the faults.

Recent Earthquakes near the San Gorgonio Pass Region

To calculate the stress interaction effects from past earthquakes, we consider the direct impact of three ground-rupturing earthquakes that occurred within the past 300 yr near the SGPr. While many smaller earthquakes have occurred within this region, these larger ground-rupturing events have the greatest potential to have impacted tractions along nearby faults.

The A.D. 1992 Landers Earthquake

The Landers earthquake occurred on 28 June 1992, rupturing five fault segments, striking northwest-southeast, in the Eastern California shear zone (Hart et al., 1993). The interaction of these faults created a linked fault network that generated an M7.3 earthquake, which is larger than expected for any single fault involved in the rupture (Aydin and Du, 1995). The total rupture length is estimated at 85 km on the primary rupture trace (Sieh et al., 1993). The epicenter was located on the southern portion of the Johnson Valley fault, and the rupture traveled northward along the Landers-Kickapoo, Homestead Valley, Emerson, and Camp Rock faults, crossing two extensional stepovers and one compressional stepover (e.g., Aydin and Du, 1995; Madden and Pollard, 2012), while only rupturing parts of the Johnson Valley, Emerson, and Camp Rock faults (Sieh et al., 1993). All of the involved faults were previously mapped, with the exception of the Landers-Kickapoo fault (Hart et al., 1993). The Johnson Valley and Landers-Kickapoo faults each slipped locally >2 m, and the central portion of the Homestead Valley fault slipped >3 m (Sieh et al., 1993; Aydin and Du, 1995). Slip exceeded 4 m on the Emerson fault, and a maximum dextral slip of ∼6 m occurred on the northern Emerson fault (Bryant, 1992, 1994; Sieh et al., 1993).

The A.D. 1812 Wrightwood Earthquake

The ∼M7.5 earthquake that occurred on 8 December 1812 (here referred to as the Wrightwood earthquake) is one of the earliest earthquakes documented in the historical records of California; the rupture origin and extent are still uncertain. Evidence of this event has been observed within several paleoseismic trench sites along the San Andreas fault north of Cajon Pass (Weldon and Sieh, 1985; Seitz et al., 1997; Biasi et al., 2002; Fumal et al., 2002; Weldon et al., 2002), with a maximum dextral slip of 4–6 m and possible northern rupture extent ∼100 km north of the Cajon Pass (Bemis et al., 2016). The southern extent of rupture is not well constrained. The most recent event recorded at Plunge Creek on the San Bernardino strand (site 3, Fig. 1) is dated to within the A.D. 1600s (McGill et al., 2002), but minor slip on secondary structures farther south along the San Bernardino strand, near Burro Flats (site 4, Fig. 1), dates to the early 1800s (Yule and Howland, 2001). Several paleoseismic sites along the northernmost strand of the San Jacinto fault record 1.8–3 m of slip during an early 1800s earthquake event (Kendrick and Fumal, 2005; Onderdonk et al., 2013, 2015). This record supports the plausibility of the 1812 earthquake jumping the <2 km extensional stepover between the San Andreas and San Jacinto faults (Fig. 1) and involving both faults (Onderdonk et al., 2013, 2015; Rockwell et al., 2015). Lozos (2016) used dynamic rupture models to investigate rupture scenarios that best fit the paleoseismic evidence and historical accounts of the Wrightwood earthquake. The models of Lozos (2016) suggest that the Wrightwood earthquake nucleated near Mystic Lake on the San Jacinto fault (site 10, Fig. 1), produced a maximum of 6 m of slip near Colton (site 9, Fig. 1), and propagated north onto the San Andreas fault (maximum of 4–5 m of slip between Cajon Pass and Wrightwood).

The A.D. 1726 Coachella Valley Earthquake

The Coachella segment of the southern San Andreas fault has not experienced a large earthquake in historical time. Paleoseismic studies reveal that the most recent earthquake is dated to A.D. 1726 ± 7 (Rockwell et al., 2018), with a possible rupture trace extending from the Salt Creek site along the Salton Sea (site 8, Fig. 1; Sieh and Williams, 1990) to the Thousand Palms site on the Mission Creek strand (site 6, Fig.1; Fumal et al., 2002), with at least 2 m of dextral offset at the Indio site on the Coachella segment (site 6, Fig. 1; Sieh, 1986).

METHODS

We use Poly3D, a quasi-static, three-dimensional boundary element method code, to simulate loading and interseismic deformation along the southern San Andreas fault system. Poly3D solves the relevant equations of continuum mechanics to calculate stresses and displacements throughout the model (e.g., Thomas, 1993; Crider and Pollard, 1998). Faults are discretized into triangular elements of constant slip (no opening or closing is permitted) within a linear-elastic half-space. The element size along the faults of the San Gorgonio Pass region (SGPr) average ∼4 km and allow for our models to capture fault irregularities as small as ∼10 km. We simulate the active fault geometry of the southern San Andreas fault, the San Jacinto fault, and the Eastern California shear zone (Fig. 2) based on the Southern California Earthquake Center’s Community Fault Model (CFM) version 4.0 (Plesch et al., 2007; Nicholson et al., 2017). The CFM is compiled from geologic mapping, seismicity, and geophysical data. While the CFM has been updated to version 5.2, the interpreted active fault geometry of the San Gorgonio Pass region is still under debate (e.g., Kendrick et al., 2015; Beyer et al., 2018; Fosdick and Blisniuk, 2018). We use version 4.0 of the CFM but include fault geometry modifications that serve both to improve the representation of the mapped active fault geometry and to improve the match of model and geologic uplift patterns and slip rates in the San Gorgonio Pass region (e.g., Cooke and Dair, 2011; Herbert and Cooke, 2012; Fattaruso et al., 2014; Beyer et al., 2018).

Faults within the CFM are defined to the base of the seismogenic crust. To simulate long-term and interseismic deformation, we extend the faults down to a freely slipping, horizontal basal crack at 35 km depth that simulates distributed deformation below the seismogenic zone (Marshall et al., 2009). Faults are likely to have distributed zones of shear rather than localized slip planes below the seismogenic crust. The deep slipping planes in the model serve as proxies for distributed loading below the seismogenic crust. This modification eliminates artifacts that develop when the long-term slip rates go to zero at the base of the CFM-defined faults (Fig. 2). Across many earthquake cycles, deformation patterns are primarily controlled by fault geometry (e.g., Dawers and Anders, 1995; Fay and Humphreys, 2005; Herbert and Cooke, 2012; Beyer et al., 2018). Beyer et al. (2018) showed that small changes in active fault configuration within the SGPr alter steady-state slip rates on local faults by as much as 6–7 mm/yr at some sites of slip rate investigation, or 35%–60% of the total slip rate. Basement rock properties would have to vary by factors of greater than three to have greater impact on the fault loading than fault geometry in this region. Therefore, to capture the first-order loading of active faults, we focus on simulating the detailed active fault configuration and do not consider potential secondary impacts of heterogeneous and/or anisotropic rock properties.

We prescribe the tectonic loading on the boundaries of the model base, far from investigated faults. Following Beyer et al. (2018), we implement an iterative technique that ensures a uniform tectonic velocity, determined from geodetic estimates (DeMets et al., 2010) at the model edges that are subparallel to the plate boundary (sides labeled I on Fig. 2) and a linear velocity gradient at the model edges that cross the plate boundary (sides labeled II on Fig. 2). The iterative approach of Beyer et al. (2018) ensures that applied velocities are within ∼1% of the desired tectonic loading at the base. Following the approach of Dorsett et al. (2019), the model uses a wide zone of applied basal velocity to ensure uniform velocity with depth at the lateral edges of the model. For faults that extend beyond our model area (San Andreas, San Jacinto, and Cucamonga–Sierra Madre fault systems), we apply slip rates to distal edge patches of these faults to prevent zero slip rates on these faults at the edge of our model. We apply 35 mm/yr dextral slip to the San Andreas fault at the northwestern edge of the model (Weldon and Sieh, 1985). At the southeastern edge, we prescribe 25 mm/yr dextral slip to the San Andreas fault and 10 mm/yr dextral slip to the San Jacinto fault (e.g., Sharp, 1981; Becker et al., 2005; Fay and Humphreys, 2005; Meade and Hager, 2005). Because of complex fault geometry and interaction among faults, deformation within the SGPr is not impacted significantly by variations in the partitioning of slip rates between the San Andreas and San Jacinto faults at this model edge (Fattaruso et al., 2014). We apply 1.6 mm/yr reverse slip (McPhillips and Scharer, 2018) to the western edge of the modeled Cucamonga fault to account for deformation along the Sierra Madre fault not included in our model.

We use a two-step modeling approach to estimate the interseismic stressing rates along the southern San Andreas fault system. The first model simulates deformation over many earthquake cycles (steady-state model) providing slip-rate information to a second model (interseismic model) that simulates the buildup of stress between earthquakes due to constant slip below the locking depth. In the steady-state model, tectonic loading is prescribed along the model edges at the base of the model, far from the investigated faults. The faults throughout the model have zero shear traction and slip freely in response to tectonic loading and fault interaction. This zero-shear traction simulates the low dynamic strength of faults during rupture (e.g., Di Toro et al., 2006; Goldsby and Tullis, 2011). We simulate interseismic deformation by applying the distribution of slip rates determined with the steady-state model to fault surfaces below the prescribed locking depth, and lock fault elements above the locking depth. The abrupt transition from locked to slipping at the specified locking depth used here produces stresses that are unreliable within one element of the transition, or ∼5 km. We use a locking depth of 25 km to ensure that our model results provide reliable fault tractions to ∼20 km depth that can be used in dynamic rupture simulations within the full depth of the seismogenic crust.

Estimating the Impact from Nearby Recent Earthquakes

We simulate the 1992 Landers earthquake, 1812 Wrightwood earthquake, and 1726 Coachella Valley earthquake by prescribing the interpreted coseismic slip distribution associated with each earthquake (e.g., Sieh, 1986; Hart et al., 1993; Onderdonk et al., 2015) to the modeled fault surfaces. We segment the rupture surface into multiple vertical segments and prescribe each segment a uniform slip according to observations at the rupture trace (Fig. 3). All other faults in the model are locked, and we do not consider the effect of tectonic loading while simulating each earthquake due to the short rupture time. The resulting static stress change due to each earthquake alters tractions along the faults within the SGPr.

Estimating Evolved Tractions

The interseismic model determines stressing rates due to deep movement below the seismogenic crust, and uses these stressing rates to calculate current shear tractions along the fault segments of the southern San Andreas fault within the SGPr using the time since the last rupture. Estimating both shear and normal tractions from stressing rates requires an assumption on how such accumulated tractions may dissipate with time. This approach relies on the premise that shear tractions that accumulate during the interseismic period are released during earthquake events (Fig. 4). Because normal tractions that accumulate in the interseismic period are not necessarily relieved upon fault slip, models of earthquake cycles require dissipation mechanisms in order to avoid singular-valued normal tractions. For example, the models of this study can estimate interseismic accumulation of normal tractions along faults within the restraining bend of the San Gorgonio Pass. However, in order to use these interseismic stressing rates to estimate present-day normal tractions, we need to know the normal traction at a specific time in the past. Furthermore, modeled coseismic slip through the San Gorgonio Pass also increases normal compression such that the accumulation of normal compression becomes unbounded within the restraining bend across multiple earthquake cycles. Estimating the evolution of normal tractions, while critical for accurate dynamic rupture models, remains a challenge for all simulations that span multiple earthquake cycles. Duan and Oglesby (2006) simulated multiple earthquake cycles by coupling a viscoelastic interseismic model with an elastic dynamic rupture model, such that normal stresses are relaxed during the interseismic period in the viscoelastic model and used as input to the dynamic rupture model. This approach can produce scenarios of potential ruptures along specific fault structures (e.g., Duan, 2019), but doesn’t simulate present-day stress because past normal stress conditions are unconstrained. Alternatively, the Rate and State Earthquake Simulator (RSQSim) employs a constant but spatially variable normal stress distribution and disregards accumulated normal tractions (e.g., Richards-Dinger and Dieterich, 2012). This approach essentially assumes that the normal stress accumulation rate everywhere equals the stress relaxation rate in the crust. Here, we follow the approach of RSQSim and do not carry the normal stressing rates through the rest of the analysis.

To estimate the shear traction that evolves over the earthquake cycle, henceforth called the evolved shear traction, we follow Smith-Konter and Sandwell (2009) and use the stressing rate information from the interseismic model and the time since last event for each fault. In this approach, we consider only large ground-rupturing events that are preserved in the paleoseismic record. The approach analyzes a coseismic stress drop corresponding to the change from static friction to dynamic friction during large ground-rupturing earthquakes (shaded region in Fig. 4). If the dynamic strength of the fault is near zero, a complete stress drop is associated with these events. Such a complete stress drop is consistent with recent field measurements of low temperatures along recently ruptured fault surfaces, a result of a very low dynamic friction (e.g., Carpenter et al., 2012; Fulton et al., 2013; Li et al., 2015), as well as high-speed laboratory frictional experiments cited above. Consequently, the associated shear traction at any time in the earthquake cycle is: 
graphic
where τ is the evolved shear traction, forumla is the shear stressing rate, and t is the time since last event. We sum the evolved shear tractions calculated by Equation 1 and the static stress changes due to nearby earthquakes to produce the present-day evolved shear tractions along the faults within the SGPr. This simplified approach to estimate the distribution of present-day shear tractions may provide more accurate initial conditions than the approach employed by dynamic rupture models of estimating tractions by resolving the remote loading onto the faults, because the evolved tractions also incorporate both loading history and fault interaction by including long-term slip rates at depth and recent nearby earthquake events (Fig. 4).

Consideration of Geometric and Tectonic Uncertainty

Using models of deformation over multiple earthquake cycles, Beyer et al. (2018) compared the slip rate distribution from six plausible active fault configuration models to available geologic slip rate data. The analysis revealed that two active fault configurations provide the best fit to the geologic observations. For this study, we use both of the two best-fit models: the inactive Mill Creek and west Mill Creek models from Beyer et al. (2018). The most pronounced fault geometry difference between the two configurations is the addition of a throughgoing Mission Creek–Mill Creek strand through the northern part of the SGPr in the west Mill Creek model. Here, we present the results of the inactive Mill Creek model geometry, and the Supplemental Material1 contains the results of the west Mill Creek model geometry. Following Herbert and Cooke (2012), Beyer et al. (2018) also tested each of the plausible fault configurations under a range of reasonable tectonic loading (plate motion of 45–50 mm/yr at 320°–325° orientation; DeMets et al., 2010); for this study, we use the mean slip rate from the end members of permissible tectonic loadings.

RESULTS

We present the interseismic stressing rates for faults of the San Gorgonio Pass region (SGPr) and show the impact of fault interaction on these rates. We analyze the results of the models that simulate three recent ground-rupturing earthquakes and the impact of these earthquakes on fault tractions within the SGPr. We then calculate the total evolved shear tractions that incorporate the impacts of both fault interaction and loading history.

Stressing Rates

Maps of interseismic shear stressing rate along the southern San Andreas fault reveal how the fault geometry controls the stressing rate distribution (Fig. 5). Figure 5 shows stressing rates for the inactive Mill Creek model configuration (Figs. 5A and 5C) and the difference in stressing rates between the two plausible active fault geometries (Figs. 5B and 5D). Dextral shear stressing rates are larger (maximum 12 kPa/yr) than the reverse-shear stressing rates (maximum ∼3 kP/yr) along the San Andreas fault. Furthermore, portions of the faults parallel to the overall plate motion, outside of the restraining bend, have a greater dextral stressing rate than faults within the bend. Dextral shear stressing rates are largest along the San Bernardino and Mission Creek strands of the San Andreas fault and decrease within the restraining bend of the SGPr (Fig. 5A). The San Gorgonio Pass thrust has an undulating strike, and small sinistral shear stressing rates occur locally along patches of the western San Gorgonio Pass thrust where the patches strike <∼265°. The reverse-shear stressing rates are near zero outside the restraining bend and increase within the bend along north-dipping fault strands that strike obliquely to the plate motion and accommodate uplift (Fig. 5C). Stressing rates increase with depth, consistent with the deep slip that is applied to faults in the interseismic model. The difference in stressing rates between the two best-fitting model geometries (Figs. 5B and 5D) are <1 kPa/yr and indicate that the west Mill Creek fault geometry produces higher dextral stressing rates (blue) throughout most of the region. The greatest difference in reverse-shear stressing rates is limited to within and just outside the bend. We only consider the inactive Mill Creek fault geometry for the rest of our analysis, and consequently, reported shear tractions may underestimate by ∼2% shear tractions if the true active fault geometry is closer to the configuration of the west Mill Creek model.

To the first order, the strike-parallel shear stressing rate along the southern San Andreas fault correlates with the orientation of the fault segments relative to the applied model loading that simulates plate motions. Previous models of the region have shown significant interaction between the San Andreas and San Jacinto faults (Herbert et al., 2014; Fattaruso et al., 2014), so we expand the analysis of stressing rates to include both of these faults in order to investigate the influence of fault interaction on stressing rates. The model produces different interseismic dextral shear stressing rates along similarly oriented portions of the San Andreas and San Jacinto faults. Figure 6 shows a gridded surface fit through model data points (white circles) of dextral stressing rates for different strikes and depths of the San Andreas (Fig. 6A) and San Jacinto faults (Fig. 6B). Stressing rates for both faults increase with depth, and in general, dextral stressing rates are higher on the San Andreas fault than on the San Jacinto fault for locations with the same strike and depth. For the San Andreas fault, maximum dextral stressing rate occurs sharply at strikes between 300°–305°. Along the San Jacinto fault, the distribution of stressing rate with fault strike shows a subdued maximum stressing rate along segments that strike ∼310°. Both of these maximum shear orientations differ from the orientation expected from resolving the regional stress tensor (black line in Fig. 6). The difference between dextral stressing rates along the San Andreas and San Jacinto faults and the expected distribution from resolved tractions demonstrates the strong impact of fault interaction on the distribution of fault stressing rates.

The impact of fault interaction is also demonstrated in the relative stressing rates on the Banning and Mission Creek strands, which differ between the two plausible fault configurations (Fig. 5). The presence of a throughgoing Mission Creek fault in the west Mill Creek model geometry (Fig. S1) shifts strike-parallel shear stressing rates from the Banning strand to the Mission Creek strand by ∼0.1 kPa/yr. These differences in stress accumulation rates over the interseismic period demonstrate that fault interaction impacts the distribution of accumulated tractions along faults within a complex system.

Impact of Stresses from Regional Earthquakes

To assess the impact of the recent nearby earthquakes along faults within the SGPr, we numerically simulate three ground-rupturing earthquakes. We simulate the Landers earthquake, Wrightwood earthquake, and Coachella Valley earthquake and examine the static stress change due to each event (Fig. 7). Because we do not consider the potential relaxation of these crustal stresses over time (e.g., Pollitz and Sacks, 2002), the fault tractions from earthquakes modeled provide an upper bound to expected tractions.

The modeled static stress change from the 1992 Landers earthquake slightly impacts tractions along faults within the San Gorgonio Pass restraining bend. The change in dextral tractions (positive) reaches a maximum of ∼0.1 MPa, along the San Gorgonio Pass thrust, and a change in sinistral tractions (negative) reaches as much as to 0.15 MPa on the southern Garnet Hill, Banning, and Mission Creek strands of the San Andreas fault (Fig. 7A). The SGPr lies in the extensional quadrant of the Landers rupture, and as a result, the fault strands within the bend are loaded with normal dip-slip tractions of ∼0.13 MPa. This dip-slip traction change effectively reduces the accumulated long-term reverse dip-slip traction on these faults.

Change in coseismic dextral tractions due to the 1812 Wrightwood earthquake increase the most (1.3 MPa) just south of the rupture limit on the San Bernardino strand, while the southernmost portion of the San Bernardino strand experiences sinistral traction changes of ∼0.25 MPa (Fig. 7B). Slip on the neighboring San Jacinto fault produces these local sinistral tractions; the southernmost portion of the San Bernardino strand did not rupture in this earthquake. The western San Gorgonio Pass thrust is loaded with dextral tractions (∼0.4 MPa), while the Garnet Hill, Banning, and Mission Creek strands experience slight (<0.1 MPa) increases in sinistral shear tractions. Furthermore, the westernmost extent of the San Gorgonio Pass thrust has normal dip-slip shear, while the rest of the thrust and the Garnet Hill strand have reverse dip-slip shear. These complex fault stressing patterns result from the location and orientation of the faults in relation to the Wrightwood rupture path. The close proximity of the dextral slip on the San Jacinto fault to the subparallel fault strands of the San Andreas fault results in sinistral coseismic traction changes on the San Andreas fault, which sits in the stress shadow of the Wrightwood earthquake.

Traction changes imposed on the San Gorgonio Pass fault strands due to the 1726 Coachella Valley earthquake reach ∼1 MPa on the Mission Creek ahead of the rupture termination and ∼0.8 on the Banning strand near the junction with the Coachella segment (Fig. 7C). Dextral tractions of up to ∼0.1 MPa extend into the restraining bend. While these nearby earthquakes are shown here to impact the SGPr, all of the resulting static stress changes due to these earthquakes are small compared to the total tractions accumulated along these faults during the interseismic period.

Estimate of Evolved Stresses

Paleoseismic data provide estimates for the time since last event, t, along active faults (e.g., Biasi and Weldon, 2009; Table 1). We estimate the total present-day traction along each fault segment by summing the tractions from Equation 1 with the static traction change of each nearby earthquake. For the San Bernardino segment, we use t from the compiled earthquake data of Biasi and Weldon (2009). Paleoseismic sites at Pitman Canyon (site 2, Fig. 1), Plunge Creek (site 3, Fig. 1), and Wrightwood (site 1, Fig. 1) provide a mean t of 207 yr for the San Bernardino segment. Paleoseismic constraints from the Thousand Palms site (site 6, Fig. 1; Fumal et al., 2002) are used for the Mission Creek strand, and from the Indio site (site 7, Fig. 1; Sieh, 1986) for the Coachella segment. These studies are in agreement that the last rupture event occurred ca. A.D. 1680. This event has been re-dated by Rockwell et al. (2018) to be ca. A.D. 1726 ± 7. For the San Gorgonio Pass thrust, Banning, and Garnet Hill strands of the San Andreas fault, we use an earthquake rupture year of ∼1400 A.D. (Heermance and Yule, 2017; Yule et al., 2014).

Due to the variable t across faults of the SGPr, the evolved shear traction distribution along the fault surfaces (Fig. 8A) differs significantly from the shear stressing rate distributions (Fig. 5). Whereas dextral-shear stressing rates are lower along the north-dipping fault surfaces within the SGPr restraining bend than on fault surfaces outside of the bend (Fig. 5A), the longer t for the faults within the restraining bend increases the total accumulated dextral shear traction within the bend relative to other faults (Fig. 8A). Similarly, although the Coachella segment and the San Bernardino strand of the San Andreas fault have greater dextral stressing rates than the restraining segment, the more recent rupture of these segments in the 1726 and 1812 events, respectively, reduces the accumulated tractions outside the bend. The largest evolved dextral shear tractions arise along the Banning and Garnet Hill strands of the San Andreas fault near the juncture with the Coachella segment of the San Andreas fault (Fig. 8A). Regions of high dextral shear traction also arise along portions of the San Gorgonio Pass thrust. The evolved reverse shear tractions are greatest along the San Gorgonio Pass thrust within the restraining bend. We note that if the true active fault geometry is more closely approximated by the alternative fault configuration (Fig. S1), the total evolved shear tractions may be underestimated by as much as 2%.

These evolved shear tractions incorporate fault interaction (Fig. 6), the t for each fault strand (Table 1), and the impact of recent nearby earthquakes (Fig. 7). To assess the impact of fault history and interaction, we compare the evolved dextral shear traction (Fig. 8A) to the fault tractions that result from resolving the regional stress tensor constrained from focal mechanism inversions onto the faults (Fig. 8B). Following Tarnowski (2017), we use the orientation of the stress field and relative magnitude of the principal stress axes from Hardebeck and Hauksson (2001) and the stress ratio, Aϕ (Simpson, 1997), of 1.5, which indicates a mixed strike-slip and thrust stress regime. The resulting stress tensor is scaled such that the change from static to dynamic friction results in a 3 MPa stress drop (e.g., Tarnowski, 2017) in order to represent the fault loading conditions preceding a large earthquake rupture. The larger magnitude of the resolved dextral shear tractions compared to the evolved tractions is due to the scaling of the regional stress to produce failure and a 3 MPa stress drop with a dynamic coefficient of friction of 0.1. The evolved tractions to the year A.D. 2019 are not explicitly at failure, and these tractions exclude those required for dynamic sliding (non-shaded region of Fig. 4). Consequently, we focus our comparison on the patterns of the evolved and resolved stresses rather than the absolute level of stress. The resolved tractions have greater lateral heterogeneity, as the tractions range from 10 MPa dextral to −5 MPa sinistral shear traction, where portions of the San Gorgonio Pass thrust receives sinistral shear from resolved loading. In contrast, the evolved tractions show dextral shear tractions everywhere on faults of the SGPr. Whereas the evolved shear tractions increase with depth, the remote stress tensor is not resolved for different depths.

DISCUSSION

Here, we discuss the impact of including fault interaction and the effects of recent nearby earthquakes on fault tractions, and the implications of this study’s findings for seismic hazard assessment.

Resolved Tractions Likely Oversimplify Initial Conditions for Rupture

Rupture propagation within dynamic models highly depends on the initial conditions used for the model (e.g., Kame et al., 2003; Duan and Oglesby, 2007; Lozos et al., 2012; Duan, 2019). These models have the power to simulate potential rupture size and extent, as well as potential rupture paths (e.g., Oglesby et al., 2003). Most rupture dynamics studies use either a homogeneous regional stress field (e.g., Lozos et al., 2012) or a regional stress field with spatially rotating principal stress (e.g., Aochi and Fukuyama, 2002) to estimate the initial shear tractions on fault segments. Resolving the regional stress tensor onto faults does not account for fault interaction or the rupture history of each fault. Whereas these effects may be minimal in regions with planar faults, we show here that within regions of fault complexity, fault interaction and loading history can advance the fault toward, or retard it away from, failure. Prescribing shear tractions that incorporate the effects of fault interaction and the loading history of each fault may improve the accuracy of dynamic rupture models.

Due to fault interactions over multiple earthquake cycles and the variable time since last earthquake event (t) across faults of the San Gorgonio Pass region (SGPr), the evolved shear traction distribution along the fault surfaces (Fig. 8A) differs from the tractions resolved from the regional stress state (Fig. 8B). Whereas resolved dextral-shear tractions are lower along the north-dipping fault surfaces within the SGPr restraining bend than on faults outside of the bend (also seen in the stressing rates in Fig. 5), the longer t for the faults within the restraining bend increases the total dextral shear traction compared to other faults (Fig. 8A). Similarly, the more recent ruptures of the A.D. 1726 and 1812 earthquake events reduce the accumulated shear tractions along faults outside of the bend, but not along faults within the bend. The largest dextral shear tractions arise along the Banning and Garnet Hill strands of the San Andreas fault, especially near the juncture of the Banning strand with the Coachella segment of the San Andreas fault (Fig. 8A). Furthermore, the evolved stresses do not produce the enigmatic left-lateral loading of the resolved stresses on portions of the San Gorgonio Pass thrust. These structures do not show surface evidence for left-lateral slip (Yule and Sieh, 2003). Consequently, the pattern of consistent dextral shear traction produced by the evolved stresses that include fault interaction and interseismic loading agrees with geologic evidence. Including evolved stress state for initial conditions within dynamic rupture models can produce more accurate assessment of rupture behavior along complex fault systems.

Within this study, we disregard the evolution of normal tractions along the San Andreas fault through the SGPr due to the challenge of producing reliable normal traction estimates. The assumption that the normal stress accumulation rate equals the relaxation rate may be reasonable along fault systems with relatively simple fault geometry, but becomes unsatisfying in regions of complex faulting. Irregularities along strike-slip faults can produce locally high accumulation rates of normal stress that outpace rates of stress relaxation. The discrepancy between rates of normal stress accumulation and relaxation accounts for many regions of nonzero uplift rate observed near active strike-slip faults, such as the San Bernardino Mountains at the San Gorgonio Pass. Because variations in normal traction along faults strongly influence dynamic rupture simulations (e.g., Aochi et al., 2002; Kame et al., 2003; Duan and Oglesby, 2007; Lozos et al., 2012; Duan, 2019), we need to better constrain these tractions to produce more accurate earthquake simulations. Resolving the details of variable normal traction along faults requires additional data sets, such as detailed stress inversions from borehole breakouts (e.g., Azzola et al., 2019), focal mechanisms (e.g., Abolfathian et al., 2019), and mechanical models of stress state that consider heat flow, topography, and rheologic variations (e.g., Luttrell and Smith-Konter, 2017). The most likely future breakthroughs will come from utilizing a combination of approaches.

Implications for Seismic Hazard

To assess how close the faults of the SGPr are to failure within our model of evolved fault tractions, we analyze the time until failure for each fault element. For this analysis, we do not explicitly consider static frictional fault strength, which would require estimates of normal traction, but instead calculate how many years are required from the present day to accumulate evolved net shear tractions equivalent to typical earthquake stress drops of 3 MPa (Fig. 9A) and 10 MPa (Fig. 9B) (e.g., Allmann and Shearer, 2009; Goebel et al., 2015). This approach estimates that fault elements on the easternmost Banning, Garnet Hill, and Mission Creek strands currently exceed 3 MPa at the base of the model (ellipse in Fig. 9A). When using 10 MPa for the accumulated traction required to trigger the next earthquake, the first fault element to fail is on the San Bernardino strand at 584 yr from now (ellipse in Fig. 9B). The difference in vulnerable position within the fault system today versus 584 yr in the future owes to the greater stressing rate along the San Bernardino strand compared to the other faults.

The estimate of time until next earthquake based on the shear stressing rates does not consider the role of fault normal tractions in earthquake initiation. Regions of higher fault-normal compression, such as within the restraining bend of the San Gorgonio Pass, may require greater accumulations of shear traction to precipitate failure than regions of lower fault-normal compression. The observed increase of earthquake stress drops with depth in a variety of tectonic settings (e.g., Trugman and Shearer, 2017; Baltay et al., 2019) confirms the normal traction dependence on earthquakes. Within the SGPr, Goebel et al. (2015) estimated much greater stress drops within the restraining bend of the San Gorgonio Pass (>20 MPa) than outside of the bend (∼4 MPa). These findings suggest that different critical levels of accumulated shear traction along different portions of the San Andreas fault system could serve as proxy for spatially variable normal traction.

The time since last ground-rupturing earthquake event for fault strands within and just outside the restraining bend are near or greater than the estimated recurrence interval for these fault strands, and paleoseismic studies in this region suggest that these faults are probably overdue, or close to failure (Philibosian et al., 2011; Rockwell et al., 2018; Onderdonk et al., 2018). Consequently, while 3 MPa represents a low stress drop in the SGPr (Goebel et al., 2015), this lower stress drop implies that these faults are close to failure, which is consistent with the earthquake clock and recurrence intervals of these faults.

CONCLUSIONS

We use three-dimensional crustal deformation models to estimate present-day fault tractions in the San Gorgonio Pass region. The models that estimate interseismic stressing rates are loaded with deep-slip rates determined from a steady-state model that simulates deformation over many earthquake cycles and explicitly includes fault interaction. Consequently, the interseismic stresses incorporate both regional tectonic loading and fault interaction. A gradient of increasing shear stressing rates with depth emerges from our models that is consistent with deep interseismic deformation. To investigate the role of fault interaction within our models, we compare our modeled dextral shear stressing rates for the San Andreas and San Jacinto faults. These faults have similar orientation, therefore the interseismic stressing rates would be similar on the two faults if based solely on orientation with respect to remote loading. Significant differences in the patterns of stressing rates along patches of the San Andreas and San Jacinto faults with similar orientations and depth arise due to the interaction between these two faults.

The total evolved present-day shear tractions along a fault include both the accumulated stressing rates since the last earthquake event and the impact of nearby earthquakes. We simulate recent nearby ground-rupturing earthquakes with coseismic models to superpose the impact of these rupture events on the stress state along the San Andreas fault within the San Gorgonio Pass region. The pattern of total evolved fault tractions differs from that of the interseismic stressing rates. Tractions are higher within the restraining bend than outside the bend because of the longer time since last event on these faults. Fault strands within the restraining bend have been loading for twice as long as the Coachella segment to the south and three times as long as the San Bernardino strand to the north. Comparison of our evolved tractions to the tractions resolved from the local stress field shows distinct differences. While the linear gradient with depth emerges from our models, the resolved tractions are not depth dependent, and the gradient must be added. Because the evolved tractions account for loading history, the largest tractions occur within the restraining bend, which has the longer time since last event.

We investigate the time needed for the accumulated net shear traction on each fault element to exceed 3 MPa and 10 MPa, typical coseismic stress drop values. Because the interseismic stressing rates differ for faults throughout the San Gorgonio Pass region, the location and timing of potential failure depends on the stress drop value used as the shear traction threshold. Assuming a lower stress drop value shows that faults in the San Gorgonio Pass are currently at failure, whereas assuming higher stress drop values does not.

This approach provides a more heterogeneous, more accurate representation of the current stress state along the southern San Andreas fault than a simple regional stress tensor. In regions of complex fault geometry such as the San Gorgonio Pass region, an “earthquake gate”, the potential for a throughgoing rupture is unclear, and stress state may have significant control on rupture behavior. Our evolved fault tractions can provide more realistic initial conditions for dynamic rupture models of these regions, and therefore improve seismic hazard assessments.

ACKNOWLEDGMENTS

The project was funded by U.S. Geological Survey (USGS) External Research Program grant G15AP00068, National Science Foundation (NSF) grant EAR 1623637, and the Southern California Earthquake Center (SCEC; contribution no. 9076). SCEC is funded by NSF cooperative agreement EAR-1600087 and USGS cooperative agreement G17AC00047. The authors thank the reviewers, Benchun Duan and Pieter-Ewald Share, for their careful reviews of this manuscript. Figure 1 was produced with the Generic Mapping Tools package (Wessel et al., 2013).

1Supplemental Material. Right-lateral and reverse dip-slip stressing rates for the alternative fault configuration presented in Beyer et al. (2018). Please visit https://doi.org/10.1130/GES02153.S1 or access the full-text article on www.gsapubs.org to view the Supplemental Material.
Science Editor: Andrea Hampel
Guest Associate Editor: Doug Yule
Gold Open Access: This paper is published under the terms of the CC-BY-NC license.