Abstract

Geologic data suggest that the Coachella Valley segment of the southern San Andreas fault (southern California, USA) is past its average recurrence time period. At its northern edge, this right-lateral fault segment branches into the Mission Creek and Banning strands of the San Andreas fault. Depending on how rupture propagates through this region, there is the possibility of a throughgoing rupture that could lead to the channeling of damaging seismic energy into the Los Angeles Basin. The fault structures and potential rupture scenarios on these two strands differ significantly, which highlights the need to determine which strand provides a more likely rupture path and the circumstances that control this rupture path. In this study, we examine the effect of different assumptions about fault geometry and initial stress pattern on the dynamic rupture process to test multiple rupture scenarios and thus investigate the most likely path(s) of a rupture that starts on the Coachella Valley segment. We consider three types of fault geometry based on the Southern California Earthquake Center Community Fault Model, and we create a three-dimensional finite-element mesh for each of them. These three meshes are then incorporated into the finite-element method code FaultMod to compute a physical model for the rupture dynamics. We use a slip-weakening friction law, and consider different assumptions of background stress, such as constant tractions and regional stress regimes with different orientations. Both the constant and regional stress distributions show that rupture from the Coachella Valley segment is more likely to branch to the Mission Creek than to the Banning fault strand. The fault connectivity at this branch system seems to have a significant impact on the likelihood of a throughgoing rupture, with potentially significant impacts for ground motion and seismic hazard both locally and in the greater Los Angeles metropolitan area.

INTRODUCTION

The Coachella Valley segment is the southernmost section of the right-lateral San Andreas fault in California (USA). Paleoseismic analysis has revealed that this segment has an average recurrence period of ∼180 yr (Philibosian et al., 2011). The most recent earthquake recorded on this section occurred ca. 1726 C.E. (Rockwell et al., 2018), and this segment is believed to represent one of the most likely San Andreas segments to rupture in the near future (Fialko, 2006; Field et al., 2015; Fumal et al., 2002; Philibosian et al., 2011). The complexity of the Coachella Valley segment increases as it approaches the eastern edge of the San Gorgonio Pass region, where the fault splits into several non-vertical and non-coplanar active strands: the Mission Creek, Banning, and Garnet Hill strands. Among these strands, the Mission Creek, which is considered to be the primary San Andreas fault strand in the southern California due to its record of large dextral displacement (e.g., Behr et al., 2010; Kendrick et al., 2015), continues northwest from the Coachella Valley toward the southeastern San Bernardino Mountains (Fig. 1). Based on geomorphologic evidence, the Mission Creek strand is the longest-lived and dominant structure accommodating relative plate motion, with a preferred slip rate of 14–17 mm/yr near the intersection with the Coachella Valley segment (Behr et al., 2010; Blisniuk et al., 2013; Matti et al., 1982; Matti and Morton, 1993; Munoz et al., 2016). However, the slip rate dies out farther west on the Mission Creek near the juncture with Mill Creek strand (Kendrick et al., 2015). On the other hand, the modern Banning strand, which continues toward the San Gorgonio Pass fault zone, has a paleoseismic slip rate of only 2.3–6.2 mm/yr, whereas the slip rate of the subparallel Garnet Hill strand is unknown (Gold et al., 2015; Heermance and Yule, 2017). However, Cardona (2016) reported that the most recent earthquake on the Garnet Hill occurred >600 yr ago.

A previous study on spontaneous rupture modeling of earthquakes initiated on the Coachella Valley segment indicates that a throughgoing rupture along the Banning strand could guide seismic energy into the Los Angeles Basin, causing significant ground motion (Olsen et al., 2006, 2008). However, Gold et al. (2015) indicated that the geologic slip rate on the Banning represents only ∼25% of the slip rate for the southern San Andreas fault (∼23–25 mm/yr), and therefore the remaining ∼75% of the slip rate could be taken up in part by other strands of the San Andreas fault—in particular, the Mission Creek strand. Thus, they inferred that a rupture originated on the Coachella Valley segment could potentially propagate along the Mission Creek strand and away from the southern California population. Understanding the most likely rupture path and the circumstances that control the rupture path is vital to determining the ability of the southern San Andreas fault system to generate large earthquakes in this region.

Recent dynamic modeling studies have elucidated the circumstances under which a rupture can propagate from one segment to another. The rupture can jump with a linking fault (Aochi and Fukuyama, 2002; Magistrale and Day, 1999; Oglesby, 2005) or without a linking fault (Harris et al., 1991; Harris and Day, 1993, 1999; Kase and Kuge, 1998, 2001; Lozos et al., 2012, 2013). Most of the studies for the case with no linking fault showed that the transfer of a rupture from one segment to another depends closely on the width and length of the stepover. For the case with a linking fault stepover, Oglesby (2005) modeled two vertical strike-slip fault segments linked with a non-vertical dip-slip fault stepover, and concluded that it is easier for the rupture to propagate to an extensional stepover with a linking normal fault than a compressional stepover with a linking thrust fault, and also that the linking fault could potentially increase the likelihood of the rupture to propagate across the stepover to generate a bigger earthquake.

This case study region of the San Gorgonio Pass hosts branched rather than stepped fault systems. The dynamics of branched faults is another case scenario of fault interaction (Aochi et al., 2000a, 2000b, 2002, 2005; Duan and Oglesby, 2005; Harris et al., 2002; Kame et al., 2003; Oglesby et al., 2003). Using dynamic rupture propagation on parallel offset faults, Harris et al. (1991) and Harris and Day (1993, 1999) showed that it is easier for the rupture to propagate to a parallel fault located in a tensional region than to one within a compressional region. However, branch fault studies have found cases where a rupture could favor the compressional branch instead of the tensional one. Aochi et al. (2000b) simulated dynamic rupture propagation on a Y-shaped branched fault system with piecewise planar segments. They observed that rupture propagation is sensitive to the initial shear stress and strongly depends on the magnitude of the angle between the branched faults; the branch fault is preferred over the continuation of the original fault for low angles. Aochi et al. (2002) also found similar results for symmetrical three-forked branched faults. They argued that when the primary fault is located in the most favorable direction (lower branching angle) with respect to a regional stress field, the rupture would prefer the compressional branch. Otherwise the tensional branch is always more favorable to rupture. Kame et al. (2003) also considered rupture on a planar fault that intersects a branch fault. They varied the angle between the main fault and the branch fault and the angle between the main fault and the maximum horizontal principal stress. They observed that for a wider angle of the principal stress (tension positive) with respect to the main fault, rupture favors the extensional branch, but as they considered shallower angles between the fault and the maximum principal stress, the compressional segment became the most favored branch for rupture. In summary, these studies showed that the ability for the rupture to jump to a branch fault and generate a bigger and more complex earthquake strongly depends on several factors, such as the applied background stress and the geometry of the different branch faults.

Here we use dynamic rupture simulations to investigate the effects of complex fault geometry on the rupture path for scenarios of earthquakes that start on the Coachella Valley segment and propagate northward. To do so, we use a finite-element model capable of simulating the dynamics of rupture propagation, slip, and seismic-wave generation based on fault geometry, fault frictional parameters, and background stress field. Our objective is to experiment with different fault geometry and pre-stress loading scenarios so that we may determine the conditions that permit the rupture to propagate (or not) along the Mission Creek, Banning, and/or Garnet Hill strands. In this way, we are able to determine which rupture scenarios are robust with respect to various unconstrained parameters, and which details of the model have controlling effects on rupture behavior.

MODEL PARAMETERS

Southern California Fault Geometry

Geophysical data do not uniquely constrain the subsurface active fault configuration of the southern San Andreas fault, so several viable interpretations have been developed (Fattaruso et al., 2014; Fialko, 2006; Fuis et al., 2012, 2017; Lin, 2013; Lin et al., 2007; Plesch et al., 2007). Consequently, we consider three alternative fault geometries (Fig. 2). For the first, we use the Southern California Earthquake Center (SCEC) Community Fault Model (CFM) version 4.0 (https://www.scec.org/research/cfm) that has been assembled using seismicity, geological mapping, and geophysical data (Nicholson et al., 2012; Plesch et al., 2007) to define the geometry of faults throughout Southern California. In this geometry (denoted model A), the Coachella Valley segment is vertical while the other segments dip to the northeast with an average angle of 80° for the Mission Creek and 70° for both the Banning and Garnet Hill strands.

Recent studies have used the fault surfaces from the SCEC CFM to simulate deformation over many earthquake cycles along the faults (Beyer et al., 2018; Cooke and Dair, 2011; Fattaruso et al., 2014; Herbert and Cooke, 2012) and compared the modeled results to Quaternary geological slip rates and uplift patterns. Despite an overall good fit of the CFM-based model to the geological data, these studies show that small discrepancies between the model and geological slip rates could be due in part to fault inaccuracy at depth where geometry is not well constrained. For example, Fattaruso et al. (2014) tested the hypothesis of Lin et al. (2007) and Lin (2013) that the Coachella Valley segment dips to the northeast. The models of Fattaruso et al. (2014) show that the dipping fault surface better reproduces the observed uplift pattern than a vertical active fault. More recently, Beyer et al. (2018) investigated alternative slip pathways throughout the San Gorgonio Pass region using a suite of six potential active fault geometries. While none of these models match the geologic strike-slip rates precisely at all the available sites, fits to the geological slip rates reveal two best-fitting models: the Inactive Mill Creek model and West Mill Creek model (Beyer et al., 2018). In this study, we consider both of these models, which are derived from modifications of the SCEC CFM. In the Inactive Mill Creek geometry, the Mission Creek is nearly vertical and the Banning, Coachella, and Garnet Hill dip to the northeast with average angles of 80°, 80°, and 60° respectively. We denote this geometry as model B. Due to the relatively shallow dip angle of the Garnet Hill strand, it intersects both the Banning and the Mission Creek, where the latter splits Garnet Hill into two separate segments: an upper and lower branch. It is based on the fault geometry model from Beyer et al. (2018).

A key difference between the West Mill Creek and Inactive Mill Creek models is the average dip angle of the Mission Creek segment. In the Inactive Mill Creek geometry (model B), the Mission Creek strand is nearly vertical, but in the West Mill Creek geometry (denoted model C), the Mission Creek has an average dip angle of 75°, so it does not intersect with the Garnet Hill strand at depth as it does for model B. Besides the variation in the dip angles, the other significant difference between models A, B, and C concerns the surface trace of those segments (Fig. 2). Both ends of the Garnet Hill strand are extended to the Banning for model A, while for models B and C, the Garnet Hill does not meet the Banning strand. There is also a clear difference in the trace of the Mission Creek for model A compared to models B and C close to the Coachella and Banning junction. In addition, at the intersection of the Mission Creek strand and Pinto Mountain fault, the Mission Creek is mapped as continuous with the Mill Creek strand (Fig. 1). Because there is no recent slip activity on the Mill Creek strand (Kendrick et al., 2015), we do not incorporate this strand in our geometry B, where the Mission Creek strand ends at the intersection with the Pinto Mountain fault. However, a recent study suggests that the westernmost part of the Mission Creek strand has had Quaternary slip (Fosdick and Blisniuk, 2018). Model geometries A and C are consistent with this interpretation: the Mission Creek extends to the Mill Creek and Galena Peak strands (Matti et al., 1983, 1985; Kendrick et al., 2015) as shown on Figure 2. In the following sections, we will use those three fault models (see Supplemental Animations S11, S22, and S33 in the supplemental material to visualize those geometries in three dimensions) to run forward models of dynamic rupture simulation and investigate the effect of fault geometry and pre-stress loading for a rupture initiated on the Coachella Valley segment and determine the conditions that permit the rupture to propagate (or not) across any of the modeled strands.

Dynamic Simulation Method

For each of the three geometries described above, we build a three-dimensional (3-D) finite-element mesh using Trelis software (www.csimsoft.com). We approximate the geometry of the fault surfaces with non-planar segments and discretize the model space with tetrahedral finite elements (Fig. 2). For all three geometries, we only allow rupture to reach a locking depth of 15 km, because local seismicity is significantly reduced below that depth (Lin et al., 2007). We choose a 200 m mesh size on the fault planes and set bias on curves extending from the fault surfaces in order to allow the mesh size to increase geometrically away from the faults. The mesh size shown on the fault segments in Figure 2 (left panel) is increased by a factor of 10 for visualization purposes. The resultant meshes for models A, B, and C have a total number of tetrahedral elements of ∼49 × 106, ∼47 × 106, and ∼42 × 106 respectively. To compute the physics involved in dynamic rupture propagation, we import these meshes into the finite-element code FaultMod (Barall, 2009), which has been tested in the SCEC–U.S. Geological Survey dynamic rupture code comparison exercise (Harris et al., 2009, 2018). In this work, we assume a semi-infinite homogeneous elastic medium where the P- and S-wave velocities are 5480 and 3160 m/s respectively and the material density is 2700 kg/m3.

Stress Field and Frictional Parameters

We test two types of background stress field: uniform (U) and regional (R) stress distribution. In our uniform model, all fault models have the same initial shear and normal stress regardless of local fault surface orientation, while in the regional model, shear and normal stress vary along the fault surfaces due to the projection of a constant regional stress field onto variable fault mesh orientations. We consider two uniform scenarios (U1 and U2), both of which are loaded uniformly with a shear stress of 10 MPa purely in the strike-slip direction. We use the strength parameter, S, to characterize these two cases. As defined by Das and Aki (1977), S is the ratio between the strength excess (yield stress minus initial shear stress) and the dynamic stress drop (initial shear stress minus sliding stress); the lower the S value, the closer the fault is to failure. Scenario U1 has an S value of 2.0, and scenario U2 has an S value of 0.3. Because the shear stress is kept fixed between the two cases, the corresponding S value is obtained by adjusting the normal stress (33.33 MPa for scenario U1 and 20.15 MPa for scenario U2). Thus, the resulting stress drop for scenarios U1 and U2 is 5 MPa and 7 MPa respectively. These two scenarios, which are chosen to demonstrate extreme cases, highlight the difference between a rupture on strong (scenario U1) versus that on weak (scenario U2) fault segments.

We also implement dynamic rupture models using the regional stress field based on the horizontal stress orientation obtained by Hardebeck and Hauksson (2001). In that study, Hardebeck and Hauksson (2001) showed that the stress field in southern California is spatially heterogeneous, with stress orientations that fluctuate between ∼N30°E to N30°W, with range of N10°E to N10°W within our study region. Therefore, we adopt two separate regional stress scenarios, in which scenario R1 has a horizontal stress orientation of N10°W and scenario R2 has a stress orientation of N10°E. We follow similar procedures as detailed in Douilly et al. (2015, 2017) to estimate the regional stress magnitude. To derive the initial shear and normal stresses on the faults for both of the stress orientations, we assume that the regional stress field is globally unchanged inside our domain, and we also apply a taper in the upper 3 km so that the shear and normal stresses decrease linearly to zero at the free surface. By doing so, we prevent a large near-surface stress drop, which is rarely observed during an earthquake. This taper is applied for both the regional stress and the homogeneous stress distributions.

Figures 3A and 3B display the resolved initial (pre-earthquake) shear and normal tractions on each fault segment from the regional stress field for geometries A, B, and C for scenarios R1 (N10°W) and R2 (N10°E). In contrast to that for the uniform stress distribution, the shear stress on the faults of these models is no longer purely strike-slip because of the non-vertical dip of the segments (see also Supplemental Fig. S24). Shear stress on the restraining bend (northwestern portion) of the Mission Creek strand is significantly reduced compared to the rest of the strand for both scenarios R1 and R2 for all three geometries (marked by the purple ellipse on Fig. 3A). This reduction is more pronounced for models A and C, where for scenario R1, shear stress on the Mission Creek strand decreases from ∼12 MPa to ∼7 MPa. The shear stress for the orientation of N10°E (scenario R2) is smaller on all segments compared to N10°W, and similar shear stress reduction (see purple ellipse location on Fig. 3A) occurs on the northwestern portion of the Mission Creek strand (∼7 MPa to ∼3 MPa). However, in contrast to the shear stress, the normal stress along that restraining bend is significantly higher (more compressive) compared to the rest of the Mission Creek, and the normal stress for the orientation of N10°E is larger on all segments compared to the N10°W scenario (Fig. 3B). This implies that it may be difficult for a rupture to break the entire Mission Creek strand in scenario R1, and difficult for the rupture to propagate on all fault segments in scenario R2.

Our frictional behavior consists of the linear slip weakening friction law (Andrews, 1976; Ida, 1972) with a Coulomb failure criteria. We assume a slip-weakening distance, static friction, and dynamic friction of 0.6 m, 0.6, and 0.15 respectively. Table 1 summarizes the values of the previously mentioned computational parameters for the dynamic rupture simulations conducted in this study.

SIMULATION RESULTS

For all of our simulations, we initiate the rupture on the Coachella Valley segment by creating a circular region of 5 km radius in which the shear stress is 10% greater than the failure stress. Figure 4 shows the distributions of final slip (40 s into simulation) for all three fault geometries for scenarios U1 and U2 (constant tractions on all fault segments). In both scenarios, the rupture propagates readily from the Coachella Valley segment to the Mission Creek strand. In scenario U1 (high strength), for all three geometries the rupture does not propagate to the Banning or Garnet Hill strands. In scenario U2 (low strength), similar to scenario U1, the rupture propagates from the Coachella Valley segment to the Mission Creek strand (see Animations S4A and S4B5). For geometry B, despite the fact that the rupture did not trigger the Banning strand, rupture did transfer to the lower branch (no slip along the upper branch) of the Garnet Hill where it intersects the Mission Creek (see red arrow in Fig. 4). A decrease in the S value from 2.0 to 0.3 causes the average slip to increase from 5 m to 7 m and also affects the rupture velocity (Animations S4A and S4B [footnote 5]). For geometry A, it takes ∼34 s for the rupture to propagate over the Coachella and Mission Creek segments for the high-strength case (S = 2.0), but for the low-strength case (S = 0.3), the rupture propagation time is ∼18 s, which is ∼2× faster. Nonetheless, the consequent increase in radiated energy is not enough to overcome the barrier for the rupture to propagate to the Banning strand. To allow the rupture to jump to the remaining strands, numerical experiments indicate that the slip-weakening distance needs to be lowered by 50% for geometries A and C and 67% for geometry B (Fig. S1 [footnote 4]). This reduction of the slip-weakening distance diminishes the number of nodes within the cohesive zone (Day et al., 2005). However, a visual inspection shows that there are still enough spatial nodes within the cohesive zone to maintain a good spatial resolution for our problem. For geometry A, the rupture propagates on both the Mission Creek and Banning segments, leaves the Banning to break the Garnet Hill entirely, then grows bilaterally on the Banning strand. For geometry B, reducing the slip weakening by 50% (i.e., 0.3 m) does allow the rupture to jump onto the Banning strand but it subsequently dies out rapidly (see Animation S56). By decreasing the slip-weakening distance, the fracture energy also decreases, increasing the radiated energy. This increase in radiated energy causes small slip to develop at the western end of the Banning strand (marked by the red circle in Fig. S1 [footnote 4]) as rupture is propagating on the Mission Creek. However, similar to the results for geometry A, the back-propagation of the rupture on the Banning segment could not break this segment entirely. Nonetheless, more slip is produced on the Mission Creek compared to the Banning and Garnet Hill strands, which implies that the Mission Creek strand is the preferred rupture path for northward-propagating ruptures under uniform loading.

Figure 5 shows the final slip distribution for scenarios R1 and R2 on all three fault geometries. Similar to the results for the constant traction models, ruptures in scenario R1 transfer easily to the Mission Creek strand, but cannot transfer to the Banning or Garnet Hill strands except on geometry B, where the rupture propagates to the lower (tensional) branch of the Garnet Hill, similar to results for scenario U2 (see Fig. S3 [footnote 4] and Animation S67). Scenario R2 also only produces rupture propagation along the Mission Creek, although it dies out shortly past the fault branch for all geometries. The near-surface small slips observed on the Garnet and Banning strands (marked by red circles in Fig. 5) are numerical artifacts due to the linear taper applied in the upper 3 km to decrease the stress to zero at the free surface, and does not have any influence on the overall rupture propagation. For both stress orientations, despite the fact that the rupture transfers easily from the Coachella Valley segment to the Mission Creek strand, it does not propagate over the compressional bend on the northwestern portion of the Mission Creek near the Pinto Mountain fault. This is due to the decrease in shear stress and increase in fault normal compression observed on Figures 3A and 3B near the compressional bend. Furthermore, the rupture dies out more rapidly for the stress orientation of N10°E compared to the stress orientation of N10°W. This is due to the fact that the shear stress magnitudes for the N10°E orientation are smaller and the normal stress higher (higher strength excess) compared to the N10°W stress orientation, which has a lower strength excess. Therefore, we could infer that the decrease in shear stress and the increase in normal stress may act as a barrier to northward-propagating rupture on the Mission Creek.

The more detailed Figures 6A, 6B, and 6C show snapshots of slip, slip rate, and shear stress change at different time steps for the stress orientation of N10°W for the respective geometries A, B, and C (see Animations S6 [footnote 7], S78, and S89). For all three geometries, the rupture propagates at sub-shear speed, and the total energy released is consistent (Mw = 7.4) among all of them. Despite the fact that the rupture does not transfer to the Banning or Garnet Hill strands, the rupture does affect the stress state on those strands. The third column of Figures 6A, 6B, and 6C shows the temporal evolution of the shear stress change (i.e., stress change with respect to the initial shear stress). As the rupture front, which is characterized by an increase in stress of ∼7 MPa, propagates from the Coachella Valley segment to the Mission Creek (time t = 8–12 s), it reduces the shear stress on the Banning and Garnet Hill by ∼2 MPa and thus pushes those two strands further from failure. This decrease in shear stress on the Banning and Garnet Hill strands is a consequence of the fact that right-lateral slip along the Coachella Valley segment relaxes shear stresses on the nearby Banning and Garnet Hill strands.

DISCUSSION

Implications of Fault Connectivity

In this study, for a northward-propagating rupture initiated on the Coachella Valley segment, we investigate the effects of fault geometry and pre-stress loading on rupture propagation and explore the most likely rupture path among several branches: the Mission Creek, Banning, and Garnet Hill strands. Our dynamic rupture simulations, for both the constant and regional stress distributions, show that the rupture is more likely to branch from the Coachella Valley segment to the Mission Creek strand instead of to the Banning strand.

An artifact of the finite-element method model may contribute to inhibiting rupture propagation beyond the branch. One challenging aspect of modeling fault interactions using most finite-element method implementations is that two or more faults cannot share common nodes to prevent a void from opening at the junction (Andrews, 1989). As an example, if we consider the triple junction in our fault system between the Coachella, Mission Creek, and Banning segments, the intersecting or edge nodes can be part of only one of those three faults. To solve this problem, instead of choosing one fault to contain the edge nodes, finite-element programs usually force slip to be zero along these nodes. Unfortunately, this option causes the fault system to be discontinuous by a gap of roughly twice the mesh element size for this particular fault system. Because the rupture initiates on the primary fault (Coachella Valley segment), it would have to jump over the barrier if the rupture propagates along any of the two branch faults (Mission Creek and Banning strands). Depending on the size of the gap, this artifact could in some cases influence the outcome of the simulation. Studies such as that of DeDontney et al. (2012) have discussed the impact of a numerical mesh gap in dynamic rupture simulation along a main planar fault with a secondary branch fault for different branching angles and S values. DeDontney et al. considered three configurations. In the first, the main fault is continuous; in the second, the branch fault is continuous; in the third, neither of them is continuous. They found that the rupture preferentially propagates on the continuous surface, and for the totally discontinuous case, the branch fault is ruptured more often than when the main fault is the continuous surface and less often than when the branch fault is continuous. Aochi et al. (2002) used a boundary-element model to test ruptures on both continuous and discontinuous fault systems. They found that for the continuous case, as the rupture front reaches the branch fault, the shear stress change dominates the fault response, while for the discontinuous case, both the shear and normal stresses could similarly affect rupture along the branch fault. Therefore, in our dynamic rupture simulations, it is possible that the gap size in this triple junction could be a contributing factor for the rupture’s lack of propagation onto the Banning and Garnet Hill strands.

To explore this issue, we modify our geometries A and C by removing the gap between the Coachella and Banning strands (referred here as models A1 and C1). Thus, we model those two segments as a single fault with a left bend, with the new Coachella-Banning segment remaining disconnected by a small gap from the Mission Creek and the Garnet Hill strands. We use those two modified geometries to investigate whether a rupture initiated on the Coachella Valley segment would still favor the Mission Creek strand instead of the Banning or Garnet Hill strand, given the greater connectivity to the Banning. We implement a dynamic rupture simulation following the previously mentioned scenarios U1, U2, and R1 on the geometries A1 and C1. Figure 7 shows the final slip distribution for scenarios U1 (S = 2.0), U2 (S = 0.3) and R1 (stress orientation N10°W) for models A1 and C1. For the constant traction model U1, the rupture continues onto the Banning strand entirely on both geometries, and triggers slip along the Mission Creek strand, but the rupture dies out rapidly. However, for scenario U2 on both geometries, the rupture propagates on both the Mission and Banning strands. On geometry A1, the rupture also triggers slip along the Garnet Hill strand but did not break this fault entirely. Despite the fact that all of the faults have the exact same initial stress conditions, one of the key differences between the A1 and C1 geometries for scenario U2 is the amount of slip on the Mission Creek strand with respect to the Banning strand. For geometry A1, more slip has been transferred to the Banning strand, while for geometry C1, the Mission Creek has more slip. This difference in slip distribution pattern is a consequence of the fault geometry. For the regional stress distribution (scenario R1), similar to results from scenario U2, the rupture also propagates on both the Mission Creek and Banning strands, but because the stress pattern on those segments depends on the orientation of the fault relative to the principal stress (in this case, N10°W), the rupture stops sooner on the Banning strand (see Animation S910 to visualize the simulations for A1 and C1). In addition, the total slip on the Mission Creek is significantly higher compared to the Banning strand for both geometries. Animation S1011 shows the strength excess for scenario R1 on geometries A and A1. The lower the strength excess magnitude, the closer the shear stress is to the failure stress. One thing to notice is that the strength excess on the Banning near the junction is overall slightly lower for geometry A1 compared to geometry A. As the rupture propagates on the Coachella Valley segment, points on the fault ahead of the rupture front experience a decrease of the strength excess to a value of ∼1.5 MPa. Even after the rupture front crosses the junction between the Coachella, Mission Creek, and Banning segments at ∼7 s, this aspect is still noticeable on the Mission Creek but not so much on the Banning segment despite the fact that the Banning and Coachella are continuous for geometry A1. In other words, even when the connectivity is changed, the overall stress favorability effect for the Mission Creek remains. Therefore, regardless of the continuity between the segments that affects the rupture propagation, the Mission Creek strand remains more favorably oriented to rupture compared to the Banning and Garnet Hill strands.

Influence of Stress Field and Geometry on Propagation Paths

The regional stress field in this study is derived from the heterogeneous stress field of Hardebeck and Hauksson (2001), where stress orientations fluctuate roughly between ∼N30°E and N30°W near our study region. Within that wide range of values, we consider stress orientations of N10°E and N10°W, and for both of those values, the rupture prefers the Mission Creek strand (the tensional branch) over the Banning or Garnet Hill strands (the compressional branches). This result is consistent with results of other studies (Kame et al., 2003; Aochi et al., 2002) that show that the tensional branch is more favorable to rupture compared to the compressional branch. However, in some cases, a rupture could favor the compressional branch. Kame et al. (2003) performed dynamic rupture modeling on two-dimensional branched fault systems with different assumptions of the pre-stress state. They inferred that the pre-stress state can have a significant impact on rupture propagation, whereby a shallow angle of the stress orientation with respect to the main fault could promote rupture propagation along a compressional branch instead of a tensional branch. To test this possibility, consistent with the relative loading of Kame et al. (2003), we apply a principal compression orientation of N40°W (scenario R3) with respect to the fault system, which results in a low angle between the stress orientation and the main fault (Coachella Valley segment), in order to investigate whether a rupture could favor the compressional Banning strand instead of the Mission Creek strand. Figure 8 shows the final slip distribution for all three geometries for the stress orientation of N40°W. For geometries B and C, the rupture indeed propagates on the (compressional) Banning and Garnet Hill strands but not the (extensional) Mission Creek (see Animation S1112 to visualize the simulation for this case). Even though the stress orientation applied to these models is outside of the ∼N30°E to N30°W stress range outlined in Hardebeck and Hauksson (2001), this result can still be considered a viable end-member scenario.

Although the strike of the Coachella Valley segment is the same for all three fault geometries, our results show that the dip angle can also influence the rupture propagation. When we initiate the rupture on the Coachella Valley segment for the stress orientation of N40°W, the rupture propagates on the compressional segments (e.g., the Banning and Garnet Hill strands) instead of the tensional segment (e.g., Mission Creek strand) for geometries B and C. However, for geometry A, the rupture dies out before reaching the junction between the Mission Creek and Banning strands. The Coachella Valley segment is vertical for fault geometry A but dips to the north for both geometries B and C, with an average dip angle of 80°. This small change in dip angle is enough to change the resolved shear and normal tractions to promote or suppress rupture propagation.

We also observe the influence on fault geometry for scenario R3 (stress orientation N40°W) where decreasing the angle between the maximum principal stress and the major fault (Coachella Valley segment) allows the rupture to propagate to the compressional segments (Banning and Garnet Hill strands) for geometries B and C (dipping Coachella Valley segment), but the rupture dies out on the Coachella Valley segment for geometry A (vertical Coachella Valley segment). Also, for geometries A1 and C1, the geometry also affects the final slip distribution (Fig. 7). Seismic hazard studies for Southern California have only considered a vertical Coachella Valley segment, but also that the preferred rupture path is through the Garnet Hill strand (Graves et al., 2008; Olsen et al., 2009). However, recent studies have produced new evidence to suggest the possibility of a northeast-dipping Coachella Valley segment (Fialko, 2006, Lin et al., 2007; Lin, 2013; Fuis et al., 2012, 2017; Fattaruso et al., 2014). As we showed in this study, a slight change in geometry (such as the dip angle of the Coachella Valley segment) could potentially affect the rupture path and thus the seismic hazard for this region. In addition, while rupture on the Banning and Garnet Hill strands instead of the Mission Creek strand is consistent with broadband ground-motion scenarios (Graves et al., 2008; Olsen et al., 2009), our scenarios imply that for a wide variety of realistic assumptions, the preferred rupture path is through the Mission Creek strand. The shorter recurrence interval of the Mission Creek strand, ∼215 yr (Fumal et al., 2002), compared to the Banning strand recurrence interval of ∼380–610 yr (Castillo et al., 2019) suggests that for any 1000 yr period, we might expect five ruptures along the Mission Creek strand and only two along the Banning strand. Consequently, both scenarios should be considered in rupture forecasting. However, ruptures along the Mission Creek for the regional case (stress orientation N10°W) do not pass beyond the region of decreased shear along the northernmost part of this strand. Thus, a rupture initiated along the southern Mission Creek might die out before reaching the northern end of this strand. This finding is consistent with geologic observations that this portion of the Mission Creek has not been active recently (Kendrick et al., 2015). However, we should emphasize that we do not take consideration of the effect of past earthquakes and stress changes due to fault interactions during the interseismic period, which could significantly influence the initial stress state (and correspond neither to a uniform traction or constant regional stress field). In addition, using earthquakes, reflection seismology, gravity, and magnetic data, Fuis et al. (2017) reported that while the shallow section of the Banning is steep, the deeper part appears to dip to the north and hence could join the Mission Creek at depths below 15 km. While none of our three geometries assume such connection at depth, this alternate configuration could allow the rupture to propagate from the Mission Creek to the Banning strand. Therefore, an important additional ground-motion scenario for Southern California could consider not only more than one configuration of fault geometry, but also a rupture that bifurcates through both the northern (Mission Creek strand) and southern (Banning and Garnet Hill) routes assuming a rupture initiation on the Coachella Valley segment.

CONCLUSION

In this study, we use dynamic rupture simulations to investigate the effects of complex fault geometry and background stress field on rupture of the San Andreas fault east of San Gorgonio Pass, with a focus on determining the conditions that could allow the rupture to propagate on the Mission Creek, Banning, or Garnet Hill strands. We consider three separate geometries based on the SCEC CFM, and assume both uniform traction and regional background stress fields. The results reveal that fault geometry has a significant impact on throughgoing rupture along this branched fault system in Southern California. Both the constant and regional stress field assumptions indicate that the rupture is more likely to branch from the Coachella Valley segment to the Mission Creek strand (a tensional branch) than the Banning or Garnet Hill strands (compressional branches). One way to obtain slip on the Banning strand would be to modify the fault geometry at depth to coincide with the model of Fuis et al. (2017) so that it joins with the Mission Creek. Our scenarios show that in a regional stress field environment, it is harder to rupture the full extent of the Mission Creek strand, and the rupture dies out rapidly for a stress orientation of N10°E compared to a stress orientation of N10°W due to the reduction of the shear stress in the northwestern portion of this segment. The suppression of northward-propagating ruptures at the northernmost portion of the Mission Creek fault may be consistent with geologic observations that this region has not had recent slip (Kendrick et al., 2015). These results have strong implications for ground motion and seismic hazard throughout heavily populated Southern California. However, we should emphasize the fact that the initial stress conditions do not take into account heterogeneity in the state of stress due to past events and fault interactions. Therefore, future studies that include stresses derived from long-term simulation could be applied to further investigate this issue.

ACKNOWLEDGMENTS

This research was supported by grants from the U.S. National Science Foundation (NSF) (award EAR-1623739), the NSF Frontiers in Earth System Dynamics program (award EAR-1135455), and a University of California President’s Postdoctoral Fellowship to RD. The figures were generated using GMT (Wessel and Smith, 1998) and ParaView (www.paraview.org). We thank Christodoulos Kyriakopoulos and Baoning Wu for fruitful discussion that helped move this project forward. We appreciated helpful comments and suggestions from Kate Scharer and an anonymous reviewer that led to improvements of this manuscript.

1Supplemental Animation S1. Animated 3-D view of geometry A. Red segment, Coachella; blue segment, Mission Creek; green segment, Banning; and orange segment, Garnet Hill. Please visit https://doi.org/10.1130/GES02192.Sa1 or access the full-text article on www.gsapubs.org to view Animation S1.
2Supplemental Animation S2. Animated 3-D view of geometry B. Red segment, Coachella; blue segment, Mission Creek; green segment, Banning; and orange segment, Garnet Hill. Please visit https://doi.org/10.1130/GES02192.Sa2 or access the full-text article on www.gsapubs.org to view Animation S2.
3Supplemental Animation S3. Animated 3-D view of geometry C. Red segment, Coachella; blue segment, Mission Creek; green segment, Banning; and orange segment, Garnet Hill. Please visit https://doi.org/10.1130/GES02192.Sa3 or access the full-text article on www.gsapubs.org to view Animation S3.
4Supplemental Figures. Figure S1: Final slip (40 s into the simulation) for fault geometries A, B, and C for the dynamic model with constant tractions on all fault segments for a slip-weakening distance of 0.3 m and an S value of 0.3—lower-strength case (scenario U2). Figure S2: Initial shear stresses for geometries A, B, and C for a dynamic rupture model from a regional stress distribution of N10°W (left column) and N10°E (right column) respectively. Figure S3: Final slip distribution on fault geometries A, B, and C for a dynamic rupture model from a regional stress distribution of N10°W and N10°E respectively. Please visit https://doi.org/10.1130/GES02192.Sf1 or access the full-text article on www.gsapubs.org to view the Supplemental Figures.
5Supplemental Animations S4A and S4B. Animation S4A: 40 s simulation of the slip distribution on geometries A, B, and C for the dynamic model with constant tractions on all fault segments for a lower-strength case (scenario U2). Animation S4B: 40 s simulation of the slip distribution on geometries A, B, and C for the dynamic model with constant tractions on all fault segments for a higher-strength case (scenario U1). Please visit https://doi.org/10.1130/GES02192.Sa4 or access the full-text article on www.gsapubs.org to view Animations S4A and S4B.
6Supplemental Animation S5. 40 s simulation of the slip distribution on geometries A, B, and C for the dynamic model with constant tractions on all fault segments for a lower-strength case (scenario 2) and for a slip-weakening distance of 0.3 m. Please visit https://doi.org/10.1130/GES02192.Sa5 or access the full-text article on www.gsapubs.org to view Animation S5.
7Supplemental Animation S6. 40 s simulation of the slip distribution on geometries A, B, and C for a dynamic rupture model from a regional stress distribution of N10°W. Please visit https://doi.org/10.1130/GES02192.Sa6 or access the full-text article on www.gsapubs.org to view Animation S6.
8Supplemental Animation S7. 40 s simulation of the slip rate on geometries A, B, and C for a dynamic rupture model from a regional stress distribution of N10°W. Please visit https://doi.org/10.1130/GES02192.Sa7 or access the full-text article on www.gsapubs.org to view Animation S7.
9Supplemental Animation S8. 40 s simulation of the shear stress change on geometries A, B, and C for a dynamic rupture model from a regional stress distribution of N10°W. Please visit https://doi.org/10.1130/GES02192.Sa8 or access the full-text article on www.gsapubs.org to view Animation S8.
10Supplemental Animation S9. 40 s simulation of the slip distribution on geometries A1 and C1 for a dynamic rupture model from a regional stress distribution of N10°W. Please visit https://doi.org/10.1130/GES02192.Sa9 or access the full-text article on www.gsapubs.org to view Animation S9.
11Supplemental Animation S10. 40 s simulation of the strength excess for geometries A and A1. All faults are represented in the first row, while the Mission Creek strand is not shown in the second row to allow for a better view of the Banning strand. Please visit https://doi.org/10.1130/GES02192.Sa10 or access the full-text article on www.gsapubs.org to view Animation S10.
12Supplemental Animation S11. 40 s simulation of the slip distribution on geometries A, B, and C for a dynamic rupture model from a regional stress distribution of N40°W. Please visit https://doi.org/10.1130/GES02192.Sa11 or access the full-text article on www.gsapubs.org to view Animation S11.
Science Editor: Andrea Hampel
Guest Associate Editor: Andrew V. Zuza
Gold Open Access: This paper is published under the terms of the CC-BY-NC license.

Supplementary data