Fault interaction is believed to influence seismicity and crustal deformation, but the mechanics of fault interaction over various time scales remain poorly understood. We present here a numerical investigation of fault coupling and interaction over multiple time scales, using the San Andreas fault and the San Jacinto fault in southern California as an example. The San Andreas fault is the Pacific–North American plate boundary, but in southern California, a significant portion of the relative plate motion is accommodated by the subparallel San Jacinto fault. We developed a three-dimensional viscoelastoplastic finite-element model to study the ways in which these two faults may have interacted (1) during and following individual earthquakes, (2) over multiple seismic cycles, and (3) during long-term steady-state fault slip. Our results show that the cluster of nine moderate-sized earthquakes (M 6–7) on the San Jacinto fault since 1899 may have lowered the Coulomb stress on the southern San Andreas fault, delaying the “Big One,” an earthquake of magnitude 7.8 or greater that may result from rupture of much of the southern San Andreas fault. In addition to the static Coulomb stress changes associated with individual earthquakes, variations of seismicity over seismic cycles on one fault can influence the loading rate on the other fault. When the San Jacinto fault experiences clusters of earthquakes such as those in the past century, the loading rate on the San Andreas fault can be lowered by as much as ∼80%. Over longer time scales, these two faults share the slip needed to accommodate the relative plate motion. Hence, an increase in slip rate on one of these two faults causes complementary decrease on the other, which is consistent with geological observations.
Fault interaction is known to play a key role in seismic activity (King et al., 1994; Harris and Simpson, 1998; Stein, 1999; Anderson et al., 2003; Lin and Stein, 2004) and evolution of geological structures (Gupta and Scholz, 2000; Peltzer et al., 2001; Li and Liu, 2007; Liu et al., 2010), yet the way in which faults mechanically couple with each other is not fully understood. Most studies of fault interaction have focused on stress transfer associated with earthquakes, through changes in the static Coulomb stress (King et al., 1994; Stein, 1999; Pollitz and Sacks, 2002). These studies show that such stress changes are largely limited to the tips of the ruptured fault segment or to the nearby faults. Long-distance fault interaction by passing seismic waves has been studied (Gomberg et al., 2001), but the durations of such dynamic stress perturbations are relatively short. Long-term mechanical coupling between faults far apart, however, has been indicated by seismic and geological evidence (Bennett et al., 2004; Dolan et al., 2007; Liu et al., 2011). The mechanisms of such fault interaction remain to be explored.
A good place to study fault interaction is southern California, where the relative motion between the Pacific and the North American plates is accommodated by the San Andreas fault and numerous subparallel faults (Fig. 1), including the San Jacinto fault, which accommodates nearly as much of the relative plate motion as the San Andreas fault does (Becker et al., 2005; Meade and Hager, 2005; Fialko, 2006), or even more (Lundgren et al., 2009). Geological slip rates on the San Jacinto fault and the southern San Andreas fault seem to be complementary over the past few million years: Increasing slip rates on the San Andreas fault were accompanied by decreasing rates on the San Jacinto fault, and vice versa (Bennett et al., 2004), indicating that these two faults are mechanically coupled.
An exploration of the ways in which these faults are coupled would contribute to a better understanding of fault interaction and earthquake hazard in southern California. The southern San Andreas fault has been seismically quiescent in the past few centuries. Its northern segments were last ruptured by the 1812 Wrightwood earthquake (M 7.5) and the 1857 Fort Tejon earthquake (M 7.9) (Fig. 1). The southernmost segments were last ruptured in 1680 by a M 7.7 earthquake (Freed et al., 2007; Jones et al., 2008) (Fig. 1). Paleoseismic data indicate ∼200–300 yr recurrence intervals on this part of the San Andreas fault (Sieh and Williams, 1990; Weldon et al., 2004), so the southern segments of the San Andreas fault are in the late phase of the seismic cycle. A rupture of much of the southern San Andreas fault could produce the “Big One,” a M ∼7.8 earthquake that would be devastating (Jones et al., 2008).
However, if the southern San Andreas fault and the San Jacinto fault are mechanically coupled, would seismicity on the San Jacinto fault affect the loading rate and stress state on the San Andreas fault? While the southern San Andreas fault has not ruptured in the past three centuries, a series of moderate-sized earthquakes (M 6–7) have frequently ruptured the San Jacinto fault since 1899 (Fig. 1; Table 1) (Bent and Helmberger, 1991; Sanders, 1993; Sanders and Magistrale, 1997). Would these earthquakes affect the “Big One” on the southern San Andreas fault? How do these two faults interact with each other over various time scales?
In this study, we use a three-dimensional (3D) viscoelastoplastic finite-element model to investigate the mechanical interaction between the San Jacinto fault and the San Andreas fault over three time scales: (1) during and following individual earthquakes on either of these faults, (2) over multiple seismic cycles, and (3) over longer terms when slip on these faults can be treated as steady state. Although some aspects of these fault interactions, especially the short-term static Coulomb stress changes associated with individual events, have been explored in various studies, mechanical coupling between faults over multiple time scales has not been investigated systematically in a self-consistent dynamic model. Using such a model, we show that, in addition to the well-studied static Coulomb stress changes associated with each earthquake, the variations of seismicity over multiple seismic cycles on one fault can significantly affect the loading rates on the other faults in the same plate-boundary fault system.
NUMERICAL MODEL OF FAULT INTERACTION
The numerical model we use here to investigate multiple time scale fault interaction is based on the viscoelastoplastic finite-element codes we have developed for simulating stress and stain evolution during seismic cycles (Luo and Liu, 2009a, 2009b) and fault interaction in a multiple fault system (Luo and Liu, 2010). The detailed viscoelastoplastic constitutive equations have been described elsewhere (Li et al., 2009; Luo and Liu, 2010). The major features pertinent to this study are described next.
To illustrate the basic mechanics of fault interaction and to avoid effects arising from irregular fault geometries, we simplify the geometry of the San Andreas fault and the San Jacinto fault as a simple bifurcation of the plate-boundary fault. Figure 2 shows the model mesh, fault geometry, and boundary conditions. The model domain includes a 20-km-thick, elastic upper crust on top of a 40-km-thick, viscoelastic lower crust and upper mantle (Fig. 2A). Using a thicker viscoelastic lower layer does not significantly affect the results (Li and Liu, 2006; Luo and Liu, 2010). Elastic modulus is 8.25 × 1010 Pa for the upper crust and 1.1 × 1011 Pa for lower crust and upper mantle; Poisson's ratio is 0.25. These are conventional values (Turcotte and Schubert, 1982); use of other reasonable values does not change the conclusions of this work. The fault zones are represented by special fault elements in a 2-km-thick, elastoplastic vertical wall (Fig. 2B) that simulates earthquakes by strain softening (see following).
The model domain is loaded by an imposed 50 mm/yr velocity on the left side, representing the Pacific–North America relative plate motion; the right side (stable North America) is fixed (Fig. 2A). On the front and back (northern and southern) sides of the model domain, motion is fixed in the x direction and free in the y and z directions (Fig. 2A). The surface is free, and the bottom of the model domain is fixed vertically but free to move horizontally.
Governing Equations and Model Rheology
We solve the stress and strain fields using the finite-element method. Our finite-element codes were developed based on the commercial software package, Finite Element Program Generator (FEPG, http://www.fegensoft.com). The finite-element codes are parallelized, and all cases were run on the University of Missouri's computer cluster (512 central processing units).
Simulating Earthquakes and Seismic Cycles
Earthquakes are usually simulated by prescribing coseismic slips in elastic or viscoelastic dislocation models (Savage, 1983; Okada, 1985; Pollitz, 2001; Freed et al., 2007), semi-analytical models (Smith and Sandwell, 2004), and using the split node method embedded in finite-element models (Dixon et al., 2002; Fay and Humphreys, 2005). These kinematical models have the advantage of being relatively simple and can be directly constrained by coseismic slips determined by geodetic and seismic data. However, the kinematical models cannot include background stress and tectonic loading, and their evolution, which are essential in our investigation of fault coupling and interaction over various time scales.
Our numerical experiments start when the model domain has been loaded to a quasi–steady state and the regional stress patterns have stabilized. At that point, the influence of the initial stress field, which is unknown and arbitrarily chosen to be zero, is minimal. When modeling specific earthquakes, we use the known seismic moment magnitude to constrain coseismic slip, which is simulated by plastic strain on the fault plane, and the amount of slip required to produce the right amount of moment. When simulating synthetic seismic cycles, we assume that an earthquake occurs when the stress reaches the Drucker-Prager yield criterion, and we use a preset stress drop to control the amount of strain softening (see following).
Earthquake-Induced Coulomb Stress Changes
Note that the change in Coulomb stress is caused by rupture of a fault plane, which perturbs the normal and shear stresses on other faults in the surrounding area. To calculate Δσcsc, one needs to know the distribution of coseismic slip, which may be inverted from geodetic or seismological data (Parsons et al., 2008; Toda et al., 2008; Shen et al., 2009). Because Δσcsc can be separated from background stress or tectonic loading, it is commonly calculated using elastic or viscoelastic dislocation models (Savage, 1983; Okada, 1985; Freed et al., 2007), which are computationally simple. For historic earthquakes for which coseismic slips are unavailable, uniform distributions of coseismic slips on faults are often assumed (Deng and Sykes, 1997a, 1997b; Freed et al., 2007).
We take a different approach here because our model needs to simulate the evolution of the full stress field for us to investigate not only the static Coulomb stress change but also fault coupling over multiple time scales. We use strain softening on the fault plane to simulate an earthquake and coseismic slip, with the plastic yield on the fault plane constrained by the moment release based on the magnitude of the event, and then we calculate the coseismic and postseismic changes in Coulomb stress (Luo and Liu, 2010). The results of this approach are comparable with those of dislocation models, validating our method (Luo and Liu, 2010). For calculating static Coulomb stress changes, our approach is no better than the dislocation models. However, with the full stress field calculated, this model allows us to simulate both the static Coulomb stress changes by individual earthquakes and the loading variations during cycles of earthquakes in a self-consistent mechanical system.
Earthquake-induced Coulomb stress changes on the southern San Andreas fault have been explored in previous studies. The last large earthquake, the 1857 M 7.9 Fort Tejon earthquake, which ruptured the northern segment of the southern San Andreas fault (Fig. 1), increased Coulomb stress of the San Jacinto fault and the southernmost San Andreas fault, and had a significant effect on seismicity patterns of these two faults (Harris and Simpson, 1996, 1998; Rydelek and Sacks, 2001; Freed et al., 2007).
Here, we focus on the Coulomb stress changes associated with earthquakes after the 1857 Fort Tejon event. Since 1899, at least nine moderate-sized (M 6–7) earthquakes have occurred on the San Jacinto fault, and only one on the southernmost San Andreas fault (Table 1; Fig. 1). Figure 3 shows the calculated Coulomb stress changes associated with the nine events on the San Jacinto fault and the one event on the San Andreas fault. The 1899 event (M 6.7) on the San Jacinto fault decreased the Coulomb stress slightly on much of the San Andreas fault, except its northern segment, where the 1948 (M 6.5) event occurred (Fig. 3A). Other earthquakes on the San Jacinto fault generally reduced the Coulomb stress on the San Andreas fault (Fig. 3B). From 1899 to the present day, the cumulative effects of the nine moderate-sized earthquakes on the San Jacinto fault were to reduce the Coulomb stress on most parts of the San Andreas fault (Fig. 3C). Freed et al. (2007) simulated large earthquakes (Mw ≥ 6.5) along the San Andreas fault system during the past two centuries and found similar results—events on the San Jacinto fault reduced the Coulomb stress on the subparallel southern San Andreas fault.
In this case, we used a viscosity value of 4 × 1019 Pa s for the lower crust and upper mantle, as suggested by previous studies of postseismic lithospheric deformation in California (Kenner and Segall, 2000; Pollitz et al., 2001). A lower value of the viscosity would predict faster decay of the postseismic Coulomb stress within this layer, and faster stress recovery within the upper crust (Li et al., 2005). The variations of viscosity values do not change the stress patterns shown here. For lithospheric deformation over longer time scales such as the cases shown here, we used a higher viscosity (∼1020 Pa s), because Pollitz et al. (2001) have shown that effective viscosity of the upper mantle could rise up to two orders of magnitude when time scales of lithospheric deformation increase, and previous mechanical models have shown that the geological slip rates on the San Andreas fault can be best fit by a viscosity of ∼2 × 1020 Pa s for the lower crust and upper mantle (Li and Liu, 2006).
Fault Interaction during Seismic Cycles
In addition to the static Coulomb stress changes associated with earthquakes, the San Andreas fault and the San Jacinto fault are expected to interact over longer time scales because they share the relative plate motion. Given that earthquakes are usually clustered and episodic, would fluctuations of seismic activity on one fault affect the loading rate on the other?
To address this question, we investigated fault interaction over multiple seismic cycles by simulating synthetic sequences of earthquakes. To better illustrate the way in which fluctuation of the seismicity on the San Jacinto fault affects loading rate on the San Andreas fault, we simulated earthquake sequences on the San Jacinto fault while not allowing the San Andreas fault to rupture. We divided the San Jacinto fault into seven segments (SJF-1 to SJF-7 in Fig. 2B), each ∼40 km long, roughly following the recent rupture history on the San Jacinto fault (Sanders, 1993; Sanders and Magistrale, 1997) and the rupture lengths of M 6–7 earthquakes (Wells and Coppersmith, 1994). When stress in one fault segment reaches the yield stress, a synthetic earthquake occurs in this segment, and the rupture is simulated with strain softening to drop 1 MPa of shear stress (Kanamori and Anderson, 1975; Smith and Sandwell, 2003) on the failed fault segment (Luo and Liu, 2010). After this earthquake terminates, and if no other fault segments reach the yield limit, the model enters an interseismic loading state until another event occurs on one or more of these fault segments. This process was repeated to simulate earthquake sequences and seismic cycles on the San Jacinto fault.
The segmental ruptures assumed here are not meant to represent real seismicity, which is more complicated (Stein and Newman, 2004; Jackson and Kagan, 2006), and inclusion of all factors that affect fault ruptures, such as discontinuities, steps, and bends, is beyond the scope of this study. Our purpose here is to see whether fluctuations of seismicity on one fault would affect the loading rate on the other.
Even with this simplified rupture pattern, the synthetic seismicity shows clusters of events separated by periods of quiescence (Fig. 4), partly because of stress triggering between events due to postseismic stress viscous relaxation (Li et al., 2007). Other factors that may contribute to earthquake clustering, such as variations of fault yield stress, damage rheology of upper crust, and healing of fault systems (Lyakhovsky et al., 2001; Kenner and Simons, 2005; DiCaprio et al., 2008), are not included here.
Our results show that the loading rate on the southern San Andreas fault varies with seismicity on the San Jacinto fault. During the quiescent periods when the San Jacinto fault stays unruptured, the loading rates on the southern San Andreas fault are ∼500–1600 Pa/yr, or ∼40%–80% higher than those of the long-term average (Fig. 5A). On the other hand, during the period when the San Jacinto fault experiences clusters of earthquakes, the loading rates on the San Andreas fault are ∼40%–80% lower than the long-term average (Fig. 5B). The spatial variations in the change of loading rates on the San Andreas fault in Figure 5 can be ascribed to two factors: The San Andreas fault in the figure includes five segments (segments C, SAF-1, SAF-2, SAF-3, and D in Fig. 2B) and is not straight, and earthquakes on the segmented San Jacinto fault have different effects on the various segments of the San Andreas fault. Where the details will vary depends on the fault segmentations in the model; these results suggest that seismic cycles on one of the faults can affect the loading rate on the other, in addition to the static Coulomb stress changes associated with each individual earthquake. In particular, the clusters of earthquakes on the San Jacinto fault since 1899 may have significantly lowered the coeval loading rate on the southern San Andreas fault.
Long-Term Fault Interaction
The San Jacinto fault and the San Andreas fault, together with other faults in the San Andreas fault system, collectively accommodate the Pacific–North America plate relative motion in southern California. Slip rates (or strain partitioning) on the San Jacinto fault and the San Andreas fault vary in different geological measurements and geodetic observations; most studies show the slip rate on the San Jacinto fault to be slightly lower than that on the San Andreas fault. Geological measurements indicate 6–23 mm/yr slip along the San Jacinto fault (e.g., Rockwell et al., 1990; Kendrick et al., 2002) and 13–28 mm/yr on the southern San Andreas fault (e.g., Weldon and Sieh, 1985; van der Woerd et al., 2006). Geodetic observations show 9–21 mm/yr slip rates along the San Jacinto fault (Fay and Humphreys, 2005; Meade and Hager, 2005; Fialko, 2006) and 22–28 mm/yr along the southern San Andreas fault (Bennett et al., 1996; Becker et al., 2005; Meade and Hager, 2005). Using combined global positioning system (GPS) and interferometric synthetic aperture radar (InSAR) data, Lundgren et al. (2009) found 24–26 mm/yr slip rate on the San Jacinto fault and 16–18 mm/yr slip rate on the southern San Andreas fault. Spinler et al. (2010) used GPS data and an elastic block model to estimate fault slip rates in this area and found ∼12 mm/yr on the San Jacinto fault and along-strike variations of slip rates on the southernmost San Andreas fault from <10 mm/yr in the north to ∼23 mm/yr in the south. Hetland and Hager (2006) argued that interseismic strain includes all effects of previous rupture events, and suggested that geodetic fault slip rates using interseismic models should be distinct from actual fault slip rates. Bennett et al. (2004) estimated slip rates on the San Jacinto fault and the parallel segments of the southern San Andreas fault in the past 5 m.y. and found complementary variations: Faster slip rates along the San Andreas fault were accompanied by slower slip rates on the San Jacinto fault, and vice versa, with the total slip rate on the two faults close to a constant of ∼35 mm/yr.
To explore long-term interaction between the San Jacinto fault and the southern San Andreas fault, we simulated synthetic seismic cycles on both faults over a longer period (10,000 yr). In this case, fault segments on both the southern San Andreas fault and the San Jacinto fault, as well as the unbifurcated fault (Fig. 2B), were allowed to rupture when the plastic yield criterion was reached. At that time, strain softening occurs on the failed fault segment to drop 4 MPa of shear stress, and the processes are repeated when other fault segments reach the plastic yield criteria. The chosen value of stress drop for each rupture was arbitrary and did not affect the main results. We changed the value from 1 MPa for simulating sequences of M 6–7 events (Fig. 4) to 4 MPa here because this experiment simulated long-term fault slip, which may include large events and aseismic creep.
The cumulative plastic slip (strain) on each fault segment yields the average long-term slip rates, and the results depend on fault properties (the plastic yield parameters). When the same fault properties are used for all fault segments, the model predicts ∼15–18 mm/yr on each of the bifurcated fault branches (the San Jacinto fault and the southern San Andreas fault) (Fig. 6A); the sum of the slip rates on these two branches is close to the ∼32 mm/yr slip rate on the single plate-boundary fault. In another case, we assumed a stronger San Jacinto fault by doubling its cohesion. This reduces the fault slip rate to ∼8–10 mm/yr on the San Jacinto fault, but increases slip rate on the southern San Andreas fault to ∼20–24 mm/yr. The sum is again close to the ∼32 mm/yr on the single plate-boundary fault (Fig. 6B).
We should emphasize that the importance here is the slip rate patterns, or strain partitioning, between these two faults and the unbifurcated plate-boundary fault. The absolute values of the predicted slip rates depend on model parameters (plastic strength of the faults), which are not well constrained. We chose the model parameters such that the predicted slip rate on the unbifurcated plate-boundary fault would be ∼32 mm/yr, which is the geological slip rate on the central San Andreas fault (http://www.consrv.ca.gov/CGS/rghm/psha/index.htm). This is lower than the rate of Pacific–North American relative plate motion (50 mm/yr in the model). In nature, this slip deficit is made up by distributed slips on other faults in the San Andreas fault system, the Eastern California shear zone (Liu et al., 2010), Baja California shear zone (Dixon et al., 2000), and some in the Basin and Range Province. Because in our model slip occurs only on the two faults, the slip deficit is expressed in elastic strain outside the faults, which is not shown here.
The results in Figure 6 are consistent with the observed complementary slip rates on the San Jacinto fault and the San Andreas fault over the past few million years (Bennett et al., 2004). Using the cumulative plastic ruptures on each segment, we also calculated moment release on each fault and found that the sum of moment release on the bifurcated San Jacinto fault and southernmost San Andreas fault is close to that on the single plate-boundary fault, consistent with the complementary moment release on the different faults of the San Andreas fault system in southern California (Dolan et al., 2007).
DISCUSSION AND CONCLUSIONS
In a complex plate-boundary zone like that in southern California, faults may interact with each other over various spatial and temporal scales. Previous models have focused on fault interaction over relatively short time scales: dynamic stress perturbations by the passage of seismic waves (e.g., Gomberg et al., 2001; Anderson et al., 2003), coseismic static Coulomb stress changes induced by earthquakes (King et al., 1994; Harris and Simpson, 1996; Stein et al., 1997; Lin and Stein, 2004), and postseismic viscous relaxation (e.g., Freed and Lin, 2001; Freed et al., 2007). Here, we have shown that over longer time scales in a mechanically coupled fault system, the fluctuation of seismicity on one fault may significantly affect the loading rate and seismicity on the other faults.
The mechanical interaction between the San Jacinto fault and the southern San Andreas fault is particularly strong because these two subparallel faults share most of the strain of the plate-boundary zone in southern California. Thus, when one fault slips more, the other has to slip less. This is consistent with the observed long-term complementary slip rates between these two faults (Bennett et al., 2004). A similar relationship has been found between the Eastern California shear zone, which accommodates up to 25% of the relative plate motion in southern California (Dixon et al., 2003; Liu et al., 2010), and the San Andreas fault system. Dolan et al. (2007) have found codependent release of seismic moment in the Los Angeles fault system and the neighboring Eastern California shear zone (Fig. 1) during the past ∼12,000 yr.
Although the model is simplified and some of the model parameters, such as the fault strength and viscosity of the lower crust and upper mantle, are not well constrained, the main model results are what one would have expected—they are the natural consequence of strain partitioning among faults in a complex plate-boundary zone. The complementary variations of Coulomb stress, loading rates, and long-term slip rates between the San Jacinto fault and the San Andreas fault all follow the principle of energy conservation: they reflect the release of strain energy, and the total energy release from all faults in a section of the plate-boundary zone has to balance with the input of strain energy, which is mainly provided by the relative plate motion.
This principle applies to fault systems in other tectonic settings, including diffuse plate-boundary zones and midcontinent regions. For example, Liu et al. (2011) have shown that the seismic moment release between major fault systems in North China is complementary, indicating that these widespread fault systems are mechanically coupled. Studying fault interaction in continental interiors is more challenging, because the tectonic loading is shared by many faults. Hence, the relatively simple fault system studied in this paper may provide useful insight into fault interaction over different time scales.
In summary, we have explored mechanical coupling between the San Jacinto fault and the southern San Andreas fault over various time scales that span from coseismic instant to multiple seismic cycles to time scales that are long enough for slip rates to be treated as steady state. Our model predicts the complementary seismic moment release and slip rates between mechanically coupled faults, as observed on the San Jacinto fault and southern San Andreas fault. Our results indicate that, in addition to the changes of static Coulomb stress between these two faults by individual earthquakes, the loading rates on one fault can be affected by seismicity on the other fault. Particularly, in addition to the decrease in the Coulomb stress on much of the southern San Andreas fault caused by the nine M 6–7 earthquakes on the San Jacinto fault since 1899, the seismic activity on the San Jacinto fault may have also significantly lowered the loading rate on the southern San Andreas fault, both contributing to delay the “Big One” earthquake in southern California.
This work is supported by National Science Foundation grants EAR-0948620 and EAR-0635574. We thank Eric Hetland, an anonymous reviewer, and Lithosphere Editor John Goodge for constructive comments and suggestions. Publication was authorized by the director of the Bureau of Economic Geology, Jackson School of Geosciences, The University of Texas at Austin.