We use three-dimensional (3-D) dynamic finite-element models to investigate potential rupture paths of earthquakes propagating along faults through the western San Gorgonio Pass, a structurally complex region along the San Andreas fault system in southern California (USA). We focus on the right-lateral San Bernardino strand of the San Andreas fault system, the oblique thrust–right-lateral San Gorgonio Pass fault zone, and a portion of the right-lateral Garnet Hill strand of the San Andreas fault system. We use the 3-D finite-element method to model rupture propagation along a fault geometry that reflects current understanding of the local geometrical complexity and is consistent with long-term loading and observed surface deformation. We test three different types of pre-stress assumptions: (1) constant tractions (assuming pure right-lateral strike-slip motion on the San Bernardino and Garnet Hill strands and oblique thrust–right-lateral strike-slip motion on the San Gorgonio Pass fault zone), (2) a uniform regional stress regime, and (3) long-term (evolved) stress from quasi-static crustal deformation modeling. Our results imply that under the more realistic regional stress and evolved stress assumptions, throughgoing rupture propagation from the southeast to northwest (i.e., from the Garnet Hill to the San Bernardino strand) may be more likely than throughgoing rupture in the reverse direction (from the San Bernardino to the Garnet Hill strand). The results may have implications for the earthquake potential in the region as well as for ground motion in the Los Angeles Basin. The results also emphasize how fault geometry and stress patterns combine to influence rupture propagation on complex fault systems.

The San Gorgonio Pass is located between the San Bernardino and San Jacinto Mountains in southern California (USA). In this region, the San Andreas fault system does not appear to have obviously continuous faults based on mapped surface traces (Fig. 1). Fault traces appear discontinuous, and many of the relationships among the various mapped traces are unclear (Allen, 1957; Matti and Morton, 1993; Matti et al., 1985; Yule and Sieh, 2003). This uncertainty in fault continuity contributes to uncertainty in the ability of the fault system to produce throughgoing ruptures and consequently large and damaging earthquakes. The entirety of the southern San Andreas fault system, especially through the San Gorgonio Pass, has not ruptured in modern seismically recorded history, but if it were to fail in a single, throughgoing rupture, large parts of southern California could experience strong shaking, including downtown Los Angeles (Olsen et al., 2006, 2008). For this reason, the San Gorgonio Pass region became a Special Fault Study Area of the Southern California Earthquake Center (SCEC), with many papers in this volume being dedicated to its earthquake potential and, in particular, the likelihood of throughgoing rupture within the San Gorgonio Pass.

Brief Geologic Considerations

The San Gorgonio Pass has a complex history of right-lateral strike-slip faulting, switching between different active fault strands, and convergence. The Mill Creek strand of the San Andreas fault system is the most geomorphically prominent fault through the pass, but its current activity is debatable (e.g., Kendrick et al., 2015; Matti et al., 1992). Many researchers have hypothesized that following the shutoff of the Mill Creek strand, the complexity of fault inter actions in the San Gorgonio Pass increased as less mechanically optimal paths to transfer slip developed (Kendrick et al., 2015; Matti et al., 1992; Yule, 2009; Yule and Sieh, 2003). However, the Mill Creek strand is commonly still incorporated into regional active fault compilations such as the SCEC Community Fault Model (Plesch et al., 2007). Modelers tend to favor simple structures in studies of hypothetical earthquake scenarios, and consequently several studies use simplified geometries that may represent older, inactive structures in the San Gorgonio Pass (Meade and Hager, 2005; Olsen et al., 2006, 2008; Smith and Sandwell, 2006). Geologic mapping and paleoseismic studies indicate that thrust faults south of the San Bernardino Mountains (Fig. 1) are likely some of the currently active structures (Allen, 1957; Kendrick et al., 2011; Matti et al., 1985; Ramzen, 2012; Scharer et al., 2013; Yule and Sieh, 2003). These faults form a less connected southern pathway for slip compared to the more connected northern pathway along the Mill Creek strand. Quasi-static deformation models have assessed alternative three-dimensional fault geometries through the San Gorgonio Pass to see which geometries matched well with geologic slip rate and uplift data (Beyer et al., 2018; Cooke and Dair, 2011; Dair and Cooke, 2009). Incorporating the Mill Creek strand reduced the match to the slip rate and uplift data compared to models with activity limited along the southern slip pathway. Consequently, we incorporate the San Bernardino strand of the San Andreas fault system, the San Gorgonio Pass fault zone, and a portion of the Garnet Hill strand of the San Andreas fault system into our semi-regional fault geometry. For brevity, we drop “of the San Andreas fault system” for the San Bernardino and Garnet Hill strands in the remainder of this paper.

Complexity in Rupture Models

Three-dimensional (3-D) dynamic rupture modeling has made a number of recent strides in incorporating realistic faulting features such as complex 3-D geometry, laboratory-based frictional parameterizations, variability in material properties, and off-fault failure. Simple models of fault stepovers (e.g., Harris et al., 1991; Harris and Day, 1993, 1999) and fault branches (e.g., Aochi et al., 2000; Kame et al., 2003; Oglesby et al., 2003) can indicate which factors, such as stepover width, branching angle, and regional stress orientation, determine the rupture path at points of geometrical complexity. Other work has investigated the intersection of strike-slip systems with dip-slip faults (Oglesby, 2005; Oglesby and Mai, 2012) and found that the nucleation location may play a key role in the likelihood of throughgoing rupture. However, models with more realistic 3-D fault geometry and heterogeneous stress may produce emergent behaviors that are rare within homogeneous and idealized models. For example, Duan and Oglesby (2007) showed that multiple earthquake cycles can greatly modify the stress near a fault branch and enable behaviors such as backward branching that seldom emerge from simple regional-stress models (Kame et al., 2003). Complex systems of faults may produce rupture behaviors that depend on the fine details of stress level and nucleation location (Anderson et al., 2003; Kyriakopoulos et al., 2019). 3-D dynamic models that incorporate more complex 3-D fault geometry, material heterogeneity, and other factors (e.g., Ando and Kaneko, 2018; Lozos et al., 2015; Ulrich et al., 2019) can produce complex behavior, although such models may rely on relatively poorly constrained assumptions about various physical parameters. In this study, we choose a moderate level of model complexity by incorporating complex 3-D fault geometry (multiple curves and branches) and multiple assumptions about fault pre-stress but with homogenous and simple parameterizations for other properties.

Previous dynamic modeling studies have investigated the ability of rupture to propagate through the San Gorgonio Pass. Olsen et al. (2008) used a relatively simple throughgoing fault structure in the pass in their dynamic models of a large earthquake on the southern segment of the San Andreas fault and found that rupture could easily propagate through the pass, channeling strong ground motion through a waveguide into the Los Angeles Basin. Douilly et al. (2020) modeled the dynamics of the branched fault structure in the eastern San Gorgonio Pass while considering rupture nucleation on the Coachella branch of the San Andreas fault system. They found that under a variety of assumptions about fault geometry and stress parameterization, rupture tended to prefer branching to the Mission Creek strand along the northern slip pathway instead of to the more southern Banning and Garnet Hill fault strands.

In all spontaneous dynamic earthquake models, the pre-stress distribution is perhaps second only to the fault geometry in controlling the rupture path and slip distribution of the system. As noted above, heterogeneous initial stresses can greatly change rupture behavior compared to simple constant traction or regional stress tensors. Long-term models using backslip assumptions (Richards-Dinger and Dieterich, 2012) and long-term loading models (Hatch et al., 2020; Stern and Cooke, 2014) can simulate the evolution of fault stresses over multiple earthquake cycles such that the fault system naturally produces distributions that are consistent with both the fault geometry and long-term loading rates. In the present work, we use the latter methodology to assign heterogeneous stresses along faults in our study region.

To perform our 3-D spontaneous dynamic rupture modeling, we use the finite-element method code FaultMod (Barall, 2009), which has been extensively verified as part of the SCEC Code Validation Project (Harris et al., 2018, 2009). It takes as its inputs the material properties, frictional formulation and parameters, fault geometry, and pre-stress field and calculates rupture propagation, stress changes, the evolution of slip and slip rate, and ground motion. For simplicity, we assume linearly elastic material behavior and homogeneous material structure (Table 1).

Mesh Generation

We utilize a fault geometry that is consistent with geologic slip rates and patterns of GPS velocities and uplift. This fault geometry is based on the SCEC Community Fault Model version 4 (Plesch et al., 2007), with small refinements introduced by Herbert and Cooke (2012), Herbert et al. (2014), and Fattaruso et al. (2014). We smooth the fault geometry from Herbert and Cooke (2012) and mesh the fault surfaces with 300 m triangular elements (Fig. 2) and the volumes with tetrahedra. A gradual increase in cell size is implemented so that elements far from the faults are three times larger than elements around the faults to minimize computational expense. While the mesh extends down to 35 km depth, we allow only the upper 20 km to be rupturable. We handle the junction line of the San Bernardino strand and portions of the San Gorgonio Pass fault zone by removing a single line of nodes along that intersection; thus, there is no slip on any faulting surface at the intersection. For simplicity, when referring to the models, we break the fault geometry down into three segments: the San Bernardino strand (SB), the western portion of the San Gorgonio Pass fault zone (WSG), and a combination of the eastern portion of the San Gorgonio Pass fault zone (ESG) and the Garnet Hill strand (GH). Breaking the faults into these sections has no geologic or physical significance; it is simply for referential ease when describing setup or results of the models. The ESG and GH sub-segments are continuous with each other. We note that our fault geometry includes a primarily vertical SB segment and north-dipping ESG and WSG segments. The GH segment makes a smooth transition from being largely vertical in the southeast to matching the dip of the ESG segment in the west, which may have implications for rupture propagation (Kyriakopoulos et al., 2017; Lozos, 2021).

Parametric and Stress-Field Assumptions

The two main factors we investigate are the effects of the different stress-field assumptions and different nucleation locations. For all models, we use a slip-weakening friction law (Andrews, 1976a; Ida, 1972; Palmer and Rice, 1973) to isolate the effects of the stress assumptions without adding frictional complexity. The results of the models are divided into three main categories based on their initial stress inputs: constant traction, regional stress, and evolved stress. In the constant traction models, the magnitude of the initial shear traction and the magnitude of the initial effective normal traction (assuming the pore pressure is included in all subsequent references to normal stress) are spatially constant on every fault element in the system. The SB and GH segments are loaded in a purely strike-slip manner, and the WSG and ESG fault portions, which are located between the SB and GH strands, have shear stress partitioned equally between right-lateral strike-slip and thrust directions (oblique pre-stress) with an instantaneous change in stress orientation at the fault-segment boundary. These simple models are denoted constant traction for comparison with our subsequent more complex and potentially realistic pre-stress assumptions. With constant frictional parameters, we assume a relative fault strength factor (Andrews, 1976b; Das and Aki, 1977),


where μstatic is the static frictional coefficient, μdynamic is the sliding frictional coefficient, σn is the initial normal stress, and τ0 is the initial shear stress, which is consistent with subshear rupture speed for a whole space (Dunham, 2007). The magnitude of the initial constant shear and normal tractions are then set assuming a 1.8 MPa stress drop (the denominator in the expression for S). We note that due to geometrically induced time-variable normal stress and dynamic overshoot, a constantinitialtraction model does not result in a constant stress drop result. Finally, after assigning stress in the above manner, we taper the shear and normal stresses linearly from their ambient value at depths below 3 km to zero at the free surface. This assumption is similar to that of Rice (1993), assuming that pore pressure increases hydrostatically down to a certain depth and lithostatically (due to disconnected pore spaces) below. We find that the homogeneous traction field in the constant traction models leads to a faster and more energetic rupture than in our other models and thus that a larger slip-weakening parameter of 0.2 m (as opposed to the 0.13 m in other parameterizations) is necessary to resolve the breakdown process. Physical and computational parameters for the constant traction models are listed in Tables 1 and 2.

For the regional stress models, we use as inputs the orientation of the tectonic stress regime and the relative magnitude of the three principal stress axes near the San Gorgonio Pass from highresolution maps of Hardebeck and Hauksson (2001). We assume a largely strike-slip–thrust stress domain with an Aϕ parameter (a shape parameter defined in Simpson, 1997, and related to the relative magnitude of the various principal stresses) equal to 1.5, as determined for depths of 10 km within the pass (Hardebeck and Hauksson, 2001). The stress tensor is oriented 10° west of north, with the principal stresses listed in Table 3, and, as with the constant traction models, we taper the regional stresses to zero above 3 km depth. The initial shear and normal stresses (resolved onto the different fault segments) are shown in Figure 3. Shear stress and normal stress are both highly heterogeneous, with the highest shear stress and lowest normal stress on the segments most aligned with the primarily strike-slip regional loading, such as the SB and GH segments, and the lowest shear stress and highest normal stress on the more thrust-directed WSG and ESG segments.

The models with evolved stress use the stress outputs of quasi-static crustal deformation modeling from Stern and Cooke (2014) and Stern (2016), based on the method used in Hatch et al. (2020). In quasi-static crustal deformation modeling, interseismic stressing rates are evaluated from plate-boundary velocities applied to the model boundaries in a two-step approach (e.g., Marshall et al., 2009). The first step finds the distribution of long-term steadystate fault slip compatible with applied plate velocities. The second step applies the steady-state fault slip below the locking depth to simulate interseismic loading on the faults above. Present-day shear-traction distributions on the various faults are estimated as the product of interseismic stressing rates and the geologic evidence for the time since the last ground-rupturing earthquake (see Hatch et al., 2020; Richards-Dinger and Dieterich, 2012; Smith-Konter and Sandwell, 2009). This simple approach assumes complete shear-stress drop during the last ground-rupturing earthquakes on the San Bernardino strand in 1812 CE (e.g., Biasi and Weldon, 2009) and the San Gorgonio Pass fault zone and Garnet Hill strand since at least 1400 CE (Scharer and Yule, 2020). This method creates a set of pre-stress conditions that should be a more realistic estimation of the pre-stress conditions along the faults in the San Gorgonio Pass region than either constant traction or a regional stress field. Hatch et al. (2020) demonstrated that in regions of complex fault interaction, such as the San Gorgonio Pass, the shear traction resulting from interseismic loading, which includes fault interactions, differs significantly from regional stresses resolved onto the same fault geometry. Stern and Cooke (2014) showed that absolute shear traction may be increased in the fault bends of the San Gorgonio Pass relative to the strands of the San Andreas fault system outside the pass. In addition to the stress interactions of active faults within the San Gorgonio Pass, the stress estimate of Stern (2016) also includes the effect of the 1992 CE Landers earthquake, which reduces shear traction along the Garnet Hill strand and increases shear traction along the San Gorgonio Pass fault zone. Due to the size difference between the elements of our long-term models and our much smaller 300 m element sizes for dynamic rupture, we interpolate our stress values to our finite-element mesh, but we do not smooth them so much as to fundamentally change them from the original output from the quasistatic modeling.

While we use the shear traction output from Stern (2016) and Stern and Cooke (2014), we do not use the normal traction output because the normal tractions are not relieved in ground-rupturing earthquakes, which complicates estimates of cumulative tractions. Obtaining the normal traction values from stressing rate would require assuming an arbitrary time for the onset of stress accumulation. Additionally, normal tractions can become unbounded near fault intersections or kinks in the elastic material. Some earthquake-cycle models implement inelastic stress relaxation or backslip, which can reduce normalstress singularities (Duan and Oglesby, 2006, 2007; Luo et al., 2020; Richards-Dinger and Dieterich, 2012). Assuming a constant normal traction or the regional normal stress in combination with the evolved shear stress would lead to many locations on the fault system being unslippable or above their yield stress at the start of the model due to the very different spatial distributions of shear and normal stress in such a model. Instead, for simplicity, we assume a constant S value of 1.32 as in the constant traction models above and use the definition of S above and the assumed frictional values to calculate our normal stress at each point on the fault. The result is a distribution of normal traction that varies spatially across the fault geometry in a manner related to the shear traction, similar to the (non-regional) stress patterns of Wen et al. (2012). Although S is held constant in the evolved stress models to obtain the initial normal-traction distribution, the stress drop is highly spatially variable on the fault system, and dynamic stress interactions during the rupture simulation enable fault strength to vary temporally at any given point. Because the stress amplitudes in the evolved stress case are a calculated result of a model, there is no additional ramp applied to the stresses in this case. The shear tractions from the model tend to decrease with distance above the locking depth (Fig. 4). The evolved stress pattern is highly heterogeneous but does not necessarily follow the regional stress pattern of Figure 3. In particular, shear and normal stress are both relatively low on the SB segment and generally higher on the WSG, ESG, and GH segments due to the differences in time since the last ground-rupturing earthquake among these faults.

It should be noted that while the shear-stress amplitudes in the evolved stress models are calculated results from a prior numerical model, the stress amplitudes in the constant traction and regional stress models are rather arbitrary and are chosen to produce slip in the range of 1–3 m, which is consistent with earthquakes having roughly the spatial extent of our models. Other than this rather rough estimation, we made no attempt to fine tune our stress patterns to give equivalent peak or average slip or seismic moment between different models. Real-world stress drops and consequent slip may differ from these values, and differences in slip amplitudes between different classes of models are not necessarily physically significant because of the radically different pre-stress patterns. Further more, earthquakes in the real world would of course not be bounded by the artificial bounds of our numerical models; thus, the final slip values should not be taken as significant in an absolute sense. The focus should be on the rupture pattern and, to a lesser extent, the spatial distribution of slip in these models.

We nucleate our preferred models by artificially increasing the shear stress to 10% higher than the yield stress (the static friction times the normal stress) in an expanding front at a speed of 0.5 times the shear-wave speed, to a critical nucleation radius. Nucleation in the regional stress and evolved stress models, due to some very lowstress regions on the fault system, requires a small slip-weakening distance (12.9 cm) and a relatively large nucleation radius (7.5 km), as displayed in Table 1; we use these values for all our preferred models, including the constant traction models. We manually inspect the time-dependent spatial pattern of stress to confirm that we resolve the slip-weakening process with at least three spatial elements within the breakdown zone. We have experimented with smaller nucleation patches of 5.0 km radius as well as models that use a decrease in static friction rather than an increase in shear traction for nucleation. In general, such models are more likely to have rupture die out shortly after the rupture leaves the forced nucleation radius, likely due to the less-energetic initial ruptures that they produce. The results of these models are shown in the Supplemental Material1. The large nucleation zone in our preferred models can be thought of as a proxy for a well-developed rupture propagating into our fault system from outside the artificial edges of the model space, which is why we choose these models for our main discussion; however, we do not claim that this nucleation is more physically plausible than other methods.

Constant Traction Models

We model three different nucleation locations in the constant traction models, one for each of the main fault segments. The resulting final slip pattern after 45 seconds is summarized in Figure 5. Apart from nucleation on the WSG, the models produce throughgoing rupture in the San Gorgonio Pass, including both directions (from south to north and from north to south) through the pass. Nucleation on the SB segment (Fig. 5A) propagates onto the ESG and GH segments but does not propagate around the acute angle to the WSG segment. Nucleation on the WSG segment (Fig. 5B) causes rupture to continue on the ESG and GH segments, but in this case, rupture does not propagate around the acute angle to the SB segment, preventing throughgoing rupture in the pass. Nucleation on the GH segment (Fig. 5C) allows rupture to propagate to the SB branch, but rupture does not propagate to the WSG; this effect is likely due to stress shadowing between the two overlapping branch segments. As will be discussed later, the constant traction results are broadly consistent with simple two-dimensional (2-D) branching fault models by Kame et al. (2003), although the physical scenario is significantly different from the earlier work. As seen in Figures S1, S4, and S7 (footnote 1), smaller nucleation patches and frictional drop–based nucleation with both large and small nucleation patches lead to results quite similar to those presented here, except that the small high-stress nucleation patch on the GH model (Fig. S1) does not cause rupture to propagate to the WSG segment.

Regional Stress Models

We find that under the assumption of a regional stress field, rupture propagation through the pass is somewhat less likely than in the constant traction case. The results for models pre-stressed using the regional stress tensor are shown in Figure 6. Rupture propagation directed from the SB strand (Fig. 6A) does not produce throughgoing rupture to the GH segment. In contrast, rupture propagation starting from the GH segment (Fig. 6C) produces throughgoing rupture to the ESG and SB segments, bypassing the WSG branch. Thus, there is a directional asymmetry for rupture propagation in this case. When the nucleation location is on the WSG segment (Fig. 6B), rupture remains confined to the same fault. We note that rupture on the ESG segment does not propagate completely to its bottom edge near the SB intersection, likely a function of the low resolved shear stress and high resolved normal stress (and thus high S) in this location (Fig. 3). The presence of many such areas unfavorable to rupture render the regional stress models more sensitive than the constant traction models to assumptions about nucleation. A smaller nucleation patch causes rupture to die out on the ESG and WSG segments prior to hitting the fault intersection, and SB nucleation does not allow propagation across the intersection (Fig. S2, see footnote 1). Nucleation based on a lowering of static friction (Figs. S5 and S8) leads to rupture dying out even earlier on the ESG and WSG fault segments and still precludes throughgoing rupture.

Evolved Stress Models

The evolved stress models shown in Figure 7 produce results somewhat different from those of either the constant traction or regional stress models. Displaying directional asymmetry like the regional stress model, the evolved stress nucleation on the GH segment produces throughgoing rupture to the SB and WSG segments (Fig. 7C), and nucleation on the SB segment (Fig. 7A) produces rupture that does not continue onto the ESG and GH segments. However, nucleation on the WSG segment (Fig. 7B) produces a small amount of throughgoing rupture propagation to the SB segment, corresponding to backward branching, which is typically not observed in simple 2-D models of branched rupture propagation (e.g., Kame et al., 2003). However, the absolute stress level, and thus the strength excess, increases as rupture propagates away from the fault intersection on the SB fault. Even with a constant S, this leads to rupture dying out before it propagates over the full extent of this segment. Nonetheless, the constant-S nature of the evolved stress models means that there is more homogeneity in rupture favorability in these models than in the regional stress case and, consequently, somewhat less dependence of the results on nucleation assumptions. Smaller nucleation (Fig. S3, see footnote 1) leads to a rupture pattern similar to that of our preferred model but with almost no rupture propagation to the SB segment in the WSG nucleation case. Nucleation via a reduction in friction and a large nucleation zone (Fig. S6) produces results similar to those of the small nucleation model above, while the reduction in friction and a small nucleation zone (Fig. S9) does not allow spontaneous rupture propagation on the SB segment.

Our results show that even on the same fault geometry, rupture path is strongly affected by assumptions about the initial stress field. Figure 8 summarizes the rupture propagation patterns for all models across all stress assumptions. We find that in two out of our three sets of initial stress assumptions, model ruptures that approach the San Gorgonio Pass from the northwest (from the SB segment) are less likely to propagate through the pass than ruptures that approach the pass from the southeast (from the GH and ESG segments). In particular, the models produce throughgoing rupture from the SB segment only in our constant traction models, which one could argue utilize the least-realistic initial stress assumptions. In both the regional stress and evolved stress cases, rupture from the SB segment dies out without propagating onto the GH segment. This lack of favorability to throughgoing strike-slip rupture from the SB segment is likely due to a combination of the static stress field being unfavorable to slip in this region and the compressional nature of the bend from the SB to the ESG and GH segments. This result is consistent with the work of Oglesby (2005), who showed that rupture around such bends between strike-slip and thrust segments may be less favorable than comparable bends between strike-slip and normal segments.

To investigate the dynamic changes in stress state, we show the change (“Delta”) in the Coulomb failure function (DCFF) on the GH-ESG and WSG segments due to slip on the SB segment (Fig. 7A) at 30.5 seconds within the evolved stress simulation (Fig. 9A). This quantity is defined as (Reasenberg and Simpson, 1992):


where Δτ is the change in shear stress from its original value and Δσn is the change in normal stress from its original value (positive in tension). Figures 9B and 9C display the separate values of Δσn and Δτ, respectively. We see that propagation in this direction corresponds to a compressional bend: slip on the SB segment induces a narrow zone of negative DCFF (blue in Fig. 9A) on the ESG segment directly to the southeast of the junction, which is due to a negative change in normal stress (clamping; Fig. 9B) that offsets a positive change of shear stress (red in Fig. 9C).

In contrast to the cases with SB-segment nucleation, earthquakes that nucleate on the GH segment, and thus approach the pass from the southeast, have the highest likelihood of propagating through the San Gorgonio Pass, with throughgoing propagation in all our models. The bend between the ESG and SB segments is extensional in this direction, aiding throughgoing propagation, which can be seen in the DCFF and Δσn calculations for the evolved stress models in Figures 9D and 9E. Of course, there are large differences in the static pre-stress levels on the different fault segments near the intersection, which also could strongly affect the directional dependence of through going rupture: rupture may propagate more easily from a slightly higher-shear-stress (and thus higher-strength-excess) region (the ESG segment) to a lower-shear-stress (and thus lower-strength-excess) region (the SB segment) than the reverse. This importance of this static pre-stress effect is further observed in that rupture propagates both directions through the pass in the constant traction models, unlike in the evolved stress models, even though both models have the same fault geometry and the same S level of 1.32.

Ruptures that nucleate on the WSG segment have the strongest sensitivity to the pre-stress assumptions but overall seem not to favor through-going propagation. Only in the evolved stress case does rupture that starts on the WSG segment connect the SB and GH-ESG segments, and in this case, rupture still dies out on the SB segment. In the regional stress case, rupture does not propagate outside of the WSG segment, and in the constant traction case, it propagates only to the GH-ESG segments, not to the SB segment. It appears that the along-strike corrugations on the WSG segment make it overall less favorable to strike-slip rupture than the somewhat smoother GH and SB segments. In the regional stress case, the more thrust-oriented portions of the WSG segment also have less resolved shear traction and more resolved normal compression, further reducing its ability to slip and propagate rupture along strike beyond the bounds of that segment.

Our results indicate that there may be a reasonably strong chance of large ruptures propagating from south to north through the San Gorgonio Pass, at least from the GH segment. Thus, our results may, to a certain extent, help to validate the assumptions of models like TeraShake (Olsen et al., 2008), which assume simpler structures in this region that facilitate throughgoing rupture from the southeast, with consequent high ground motion in the Los Angeles Basin due to waveguide and directivity effects. It is important to note that while we believe that evolved stresses have great promise in more accurately modeling the rupture dynamics of geometrically complex fault systems, they do introduce some uncertainty, particularly in the implementation of normal traction in the current methodology. Of course, all assumptions of fault traction have large amounts of uncertainty, so ideally, multiple realizations of pre-stress field should be considered before drawing conclusions.

It should be noted that our models only investigate the western portion of the San Gorgonio Pass, with nucleation on the GH segment acting rather like a proxy for rupture being fed in from the portion of the San Andreas fault system to the southeast of our model region (e.g., the Coachella segment). Our work does not address the likelihood for whether rupture from the Coachella segment of the San Andreas fault system would feed into our model space via the GH segment. Measurements of slip rate on nearby fault segments (Behr et al., 2010; Blisniuk et al., 2013; Matti et al., 1982; Matti and Morton, 1993; Muñoz Zapata, 2017) imply that rupture would be more likely to propagate to the Mission Creek rather than the Banning and Garnet Hill strands at their intersections east of the San Gorgonio Pass, which would mean that there is not a high likelihood of rupture being fed into our model space from the southeast. The 3-D dynamic models of Douilly et al. (2020) also imply similar conclusions. Therefore, even though many of our models produce throughgoing rupture from the southeast to the northwest through the San Gorgonio Pass, this does not necessarily mean a high likelihood of a TeraShake-style rupture propagating over most of the southern San Andreas fault system. Rupture nucleating on the GH segment in our models may still propagate through the pass, but directivity toward the Los Angeles Basin would be reduced, with potential implications for ground motion therein.

Although our faulting configuration is considerably more complicated than the 2-D branching models of Kame et al. (2003), it is instructive to compare our results to this prior modeling work. Our models produce more complex behaviors than the 2-D models of Kame et al. (2003), such as potential backward branching and rupture termination at a simple bend, which could be facilitated by the small separation of the fault segments at the junction (Fliss et al., 2005), but are also likely due to the complex stress field in the vicinity of the branch. Multi-cycle modeling by Duan and Oglesby (2007) show how a variety of rupture behaviors at 2-D fault branches can be facilitated by accounting for the long-term earthquake history, and our models (particularly the evolved stress cases) share these characteristics. Of course, with the capacity of rupture to propagate in 3-D, and with highly variable shear pre-stress direction, it is difficult to draw broad-scale conclusions about faulting behavior. Our results show the value of running multiple models with multiple stress assumptions to determine a plausible range of possible faulting outcomes in future earthquakes.

It is important to note that our dynamic modeling results for potential earthquake rupture propagation in the western San Gorgonio Pass are not comprehensive. Most importantly, there is no true consensus on the fault geometry in this region, and we model earthquakes with a single active-fault geometry. Fault connectivity could greatly impact rupture propagation, and experiments with different assumptions about fault geometry (such as in Douilly et al. [2020]) are clearly necessary before drawing firm conclusions. The importance of assumptions about the frictional parameterization and nucleation should also not be underestimated. As shown in the Supplemental Material (footnote 1), different assumptions about the radius of the artificial nucleation patch size and method of nucleation can affect the ability of rupture to propagate in the region. In our preferred models, we use a 7.5 km nucleation radius as a proxy for a more fully developed rupture propagating inward from outside the artificial bounds of our model space; with smaller nucleation sizes, rupture is much less likely to propagate across the San Gorgonio Pass in either direction, with the exception of the constant traction models. Using a nucleation method via a reduction in static friction rather than an increase in shear stress in the nucleation zone similarly reduces the ability of rupture to propagate within and between fault segments. Informal experiments with different frictional coefficients and slip-weakening distances (not shown in this current work) can result in different rupture-propagation patterns, with a general observation that higher slip-weakening distance hinders rupture propagation at fault discontinuities, consistent with work on fault stepovers by Lozos et al. (2014). None of these alternative results should be discounted, given the uncertainties in input parameters.

Taken together, our results imply that the ability of rupture to propagate through the San Gorgonio Pass depends very strongly on the pre-stress assumptions and rupture propagation direction (i.e., nucleation location). We have modeled the dynamics of potential earthquakes under assumptions of constant traction, a uniform regional stress field, and a stress field drawn from a long-term model of tectonic loading, fault slip, and past recorded earthquakes. In general, it appears that across multiple assumptions about fault pre-stress, throughgoing rupture propagation from the southeast to northwest may be somewhat more likely than throughgoing rupture in the reverse direction, although that statement implicitly assumes that rupture has an equal chance of arriving from the direction of the Garnet Hill and Banning strands as from the San Bernardino strand, which may not be the case (e.g., Blisniuk et al., 2013). Additional models from the ones explored here, with alternative fault geometries, will be necessary to draw firmer conclusions and statistical significance. Nonetheless, the results may have implications for potential ground motion in the Los Angeles Basin, which may experience greatly amplified ground motion from throughgoing south-to-north rupture due to waveguide and resonance effects (Olsen et al., 2006, 2008). Our results imply that such rupture is possible, although more work is required to determine the likelihood of rupture propagating across the multiple geometrical discontinuities between the Coachella segment and the San Bernardino segment, as well as to determine the implications for ground motion in Los Angeles from ruptures initiating on the Garnet Hill segment rather than farther to the southeast.

This work also shows that it is possible to incorporate long-term evolved stresses into dynamic rupture modeling. Although the evolved stress model is one of three assumptions we make in the current paper, we believe that this form of long-term modeling to obtain reasonable pre-stress patterns for dynamic modeling has great potential and should be incorporated into dynamic models in which real-world rupture propagation potential is being explored.

The authors gratefully acknowledge John Matti and Katherine Kendrick for their help and perspective in this work and for producing this paper’s Figure 1 based on their earlier work. We also are grateful for helpful conversations with Doug Yule, Kate Scharer, and numerous SCEC researchers working on the San Gorgonio Pass Special Fault Study Area. The paper was greatly improved under review by suggestions from Alex Hatem and an anonymous reviewer. This work was supported in part by U.S. National Science Foundation grant EAR-1623739 to Oglesby and Cooke, U.S. Geological Survey–National Earthquake Hazards Reduction Program grant G15AP00068 to Oglesby and Cooke, and SCEC grant 13130 to Oglesby and Cooke. This research was supported by the Southern California Earthquake Center (Contribution No. 11873). SCEC is funded by NSF Cooperative Agreement EAR-1600087 and USGS Cooperative Agreement G17AC00047.

1Supplemental Material. Contains alternate dynamic modeling results using different nucleation methods. Please visit https://doi.org/10.1130/GEOS.S.20493282 to access the supplemental material, and contact editing@geosociety.org with any questions.
Science Editor: Andrea Hampel
Associate Editor: Craig H. Jones
Gold Open Access: This paper is published under the terms of the CC-BY-NC license.