“Aftershock Faults” and What They Could Mean for Seismic Hazard Assessment

Systems of faults collectively accommodate plate boundary strain, because it is rare that any one fault can be perfectly aligned with relative plate motions (e.g., Fitch, 1972; Michael, 1990; Gauriau and Dolan, 2021; Okuwaki and Fan, 2022). Terms such as strain partitioning and/or accommodation zones describe a process in which slip on a fault misaligned to some degree with relative plate motions causes elastic strain accumulation in the surrounding crust that must ultimately break new faults to be relieved. At smaller scale, fault bends and stepovers require similar accommodation. In this article, we examine dependence amongst faults in the California system. We calculate the degree to which faults are stressed directly by relative plate motions versus those that receive no plate‐motion stress, and we also calculate interaction stresses on faults that result from slip on the dominant fault of the California system—the San Andreas.

Our hypothesis is that some accommodating structures are only activated by earthquake ruptures on other faults, and otherwise lie fallow. We know from global observations (e.g., Omori, 1894; Dieterich, 1994; Parsons, 2002) that a sudden stress increase on neighboring faults caused by a mainshock leads to aftershocks that obey an Omori law rate decay with time. If the stressing process of a fault only occurs by interactions with other faults, it can then be thought of as an “aftershock fault.” This classification could become important in time‐dependent earthquake rupture and hazard forecasts (e.g., Field et al., 2014, 2015), because aftershock faults are not independent seismic sources. This article is not an exhaustive analysis but rather is intended to illustrate a spectrum of fault loading conditions, and how they might relate to earthquake rupture forecasting in California and generally.

We make estimates of tectonic and interaction loading on all faults within the California system, as identified by Dawson and Weldon (2013) for use in earthquake rupture forecasting to explore the relative influence of these two stressing mechanisms. We begin by examining stressing from relative motions between the Pacific plate and Sierra Nevada block (PAC‐SN). We choose this reference frame over the broader Pacific–North American motions, because the Sierra Nevada block moves relative to North America and has little internal deformation (e.g., Dixon et al., 2000). Considerable Pacific–North American plate shear occurs east of the Sierra Nevada (e.g., Parsons and Thatcher, 2011), meaning that Pacific–North American relative motion directions are a less accurate characterization of strain west of the Sierra Nevada block.

We can then make a simple assessment of fault loading by comparing the alignment of PAC‐SN relative plate motions with fault slip directions. Parsons et al. (2013) computed fault hanging‐wall displacement directions and rates for every California fault subsection based on consensus values of strike, dip, and rake (Dawson and Weldon, 2013) that are projected into the horizontal plane (Fig. 1a). We make a comparison of the unit vectors of fault slip direction and relative plate‐motion direction by calculating their dot product, in which a resulting value of 1 implies perfect alignment, a value of 0 implies 90° misalignment, and −1 implies a 180° difference (Fig. 1c). A mapping of the dot products shows that most California faults are well aligned with relative plate motions, though there are some notable exceptions including the Garlock fault, some Transverse Range faults, and the Great Valley fault system (Fig. 1c).

We make these calculations for all California faults with the understanding that our PAC‐SN block‐loading model does not apply to every fault; for example, results for faults east of the Sierra Nevada may have little meaning, because these faults would most likely be loaded by relative North America–Sierra Nevada plate motions. We do not want to use a subjective, edited selection of faults, so we include them all. The use of hanging‐wall displacement directions is an arbitrary choice, as footwall displacements could have been chosen instead. Displacement directions depend on fault dip directions, strikes, and rakes; the choice of hanging wall versus footwall is only an issue for some east–west‐trending dip‐slip faults that have displacement directions that are nearly parallel to relative motion vectors but can have negative dot products, because the hanging‐wall directions are ∼180° off relative motion vectors. This happens with some thrust faults in the Transverse Ranges (Fig. 1c), which have been individually corrected. The stressing rate calculations discussed next do not have this issue.

We make an additional independent estimate of fault loading by calculating the stress accumulation on California fault planes from relative PAC‐SN plate motions. We do this with a finite‐element model of the California elastic crust that is 20 km thick and made up of 20 node 10 km cubic elements (Fig. 2a). We opt to conduct continuum models for stress accumulation on faults instead of a back‐slip model (Savage, 1983), in which buried dislocations with opposite rakes are slipped beneath faults to load them. The back‐slip model assumes the existence of deep loading sources that are aligned directly with faults, even if those fault displacements are misaligned with relative plate motions. Because our main goal is to reconcile fault orientations with relative plate motions, we conclude that the best approach is to model the stress field as resulting from elastic strain in the crustal volume from relative plate motions.

Therefore, the stress tensor resolved onto the plane is given by $σs′=Aσs$⁠, in which the components of interest from this matrix are those that define the traction (T) on the fault plane as $T=[s13,s23,s33]$⁠. The along‐strike (⁠$s13$⁠) and downdip (⁠$s23$⁠) shear components resolved on the fault surfaces are compared with the fault subsection rake components and are net positive if the stress vector sums have any positive value in the rake directions. Otherwise, if there is no positive rake‐parallel component of the shear stress, then the negative value in the rake direction is reported (Fig. 2b). We interpret zero or negative tectonic stressing as an indication that either these fault subsections have incorrect rakes assigned or that they are not directly stressed by relative plate motions. An example test shows that the modeled stress‐loading profile on the San Andreas fault has a similar trend to the observed slip rate profile (Fig. S1, available in the supplemental material to this article).

Modeled tectonic stressing results show similarities with the simple dot product test (Fig. 1c), with most faults within the broader San Andreas system receiving positive tectonic stress loading (Fig. 2b). By “San Andreas system,” we mean the set of approximately parallel major right‐lateral strike slip faults that are well aligned with relative PAC‐SN motions. Examples of large, poorly aligned faults that receive no tectonic stressing according to the model include faults within the Great Valley thrust system, the Garlock fault, and east–west‐trending Transverse Range oblique strike‐slip faults (Fig. 2b). The existence of active faults that are apparently not loaded by tectonic motions implies that they are stressed by other means. We thus investigate secondary loading by fault interactions as a potential stressing mechanism for these faults.

We calculate stress transfer to all California faults from right‐lateral strike slip on the San Andreas fault, which is the dominant fault within the California system. The San Andreas fault is mostly well aligned to relative plate motions (Fig. 1), has the highest long‐term slip rate (e.g., Dawson and Weldon, 2013), and receives the largest stress loading from relative plate motions (Fig. 2). There are multiple parallel faults in the San Andreas system that also act to accommodate relative plate motions; slip on these faults has a similar stress transfer effect onto subfaults because of their similarity in strike and rake to the San Andreas system. Thus, stress transfer from San Andreas fault slip can be thought of as a generalization for slip within the broader right‐lateral system.

We create elastic dislocations from the UCERF3 fault subsection database (Dawson and Weldon, 2013), with all faults except the San Andreas treated as receiver faults that have defined strikes, dips, and rakes taken from geologic estimates. The San Andreas fault subsections are slipped in proportion to their relative long‐term slip rates in a boundary element model (Okada, 1992), and shear stress changes are resolved on all other California faults relative to their strikes, dips, and rakes (Fig. 3a). We calculate shear stress changes rather than Coulomb stresses, because shear stresses drive slip and are most comparable to our calculation of tectonic stresses. In addition, the normal stress term in Coulomb calculations is modulated by fault friction coefficients, which are subject to considerable uncertainty.

A closer look (Fig. 3a,b) reveals that many of the faults that are not loaded by modeled relative plate motions because of their orientations (Figs. 1 and 2) are stressed by slip on the San Andreas fault. These include the Great Valley thrusts, the Garlock fault, and multiple faults within the Transverse Ranges and Coast Ranges (Fig. 3). In general, it appears that these faults accommodate stresses resulting from slip through significant bends on the San Andreas fault, such as the westward‐directed “Big Bend” in southern California, and on San Francisco Peninsula. The Garlock fault might also be loaded by secondary interactions such as relative motion between the Mojave and Sierra Nevada blocks. In addition, there is a slight westward trend of the central San Andreas fault (the Carrizo, Cholame, and Creeping segments) of about 12° compared with relative PAC‐SN plate motions that produces compression and loads faults in the Great Valley thrust system.

We see that many of the faults that are misaligned with relative plate motions (Fig. 1), and/or that are not stressed by relative plate motions (Fig. 2), are also those that have positive stress loading from interactions with the San Andreas fault (Fig. 3). We interpret these faults as “aftershock faults,” meaning that they most likely slip during periods after earthquakes that have occurred on the San Andreas or San Andreas parallel faults. It has been observed globally (Parsons, 2002) that earthquake rates on faults after positive stress interactions with mainshocks follow Omori law (⁠$time−1$⁠) decays (e.g., Omori, 1894; Dieterich, 1994). Faults that fall into this category are mapped in Figure 4a.

We separate the spectrum of fault‐loading modes into four categories (Fig. 4c): (1) faults that are loaded primarily by relative PAC‐SN motions and are not loaded by interactions with the San Andreas fault, (2) faults loaded by interactions with the San Andreas fault and relative PAC‐SN motions, (3) faults that are loaded primarily by interactions with the San Andreas fault system and are not loaded by relative PAC‐SN motions, and (4) faults that are not loaded by interactions with the San Andreas fault nor relative PAC‐SN motions. We assess the relative importance of these by the number of fault subsections (n in Fig. 4c), and their relative cumulative annual moment release calculated by reported long‐term subsection slip rates and areas (Dawson and Weldon, 2013). The total annual moment release is restricted here to mapped faults.

Most California fault subsections (1611/2614) and annual moment (74%) fall within the first category that are primarily loaded by relative PAC‐SN motions. These tend to be major right‐lateral strike‐slip faults like the San Andreas and its many parallel strands. This is also illustrated in Figure 4d, which shows that most faults that are well aligned with relative plate motions (dot products ∼1) have negative stress change values. The next highest category in terms of subsections (413) and annual moment release (18%) are faults loaded both by interactions with the San Andreas fault and relative PAC‐SN motions. The third category consists of aftershock faults, which are only loaded by interactions with the San Andreas fault; this group consists of 337 subsections and contributes 4% of the annual moment release (these faults are mapped in Fig. 4a). Finally, the group of faults that is loaded neither by interactions with the San Andreas fault nor relative PAC‐SN motions consists of 253 subsections and releases 4% of the annual moment. These faults are mapped in Figure 4b, and it is apparent that most of them are located outside of the zone of relative PAC‐SN motion loading, lying north of the Mendocino triple junction, and are thus stressed by Cascadia subduction zone convergence. In addition, many of the other faults in this category are located east of the Sierra Nevada block, meaning that a PAC‐SN relative motion‐loading model does not apply to them. There are a few (Fig. 4b) that lie within the realm of PAC‐SN relative motion loading. These faults must be loaded by some other mechanism, and their presence suggests uncertainties with our analysis, either in assigned fault geometries, relative plate‐motion directions, or interaction loading from faults other than the San Andreas.

There were six historical M≥6 earthquakes that occurred near aftershock faults (Table S1). We investigate each to assess consistency with our model (Fig. S2). We find that three of these events were triggered by other earthquakes, two were actually located on or near tectonically loaded faults, and one was inside the Sierra Nevada block as modeled (see the supplemental material for more details).

There are multiple approaches to seismic hazard assessment. Standard probabilistic seismic hazard assessment (PSHA) treats seismic sources as following a Poisson process (random in time). These sources are convolved with ground‐motion prediction equations to estimate the probability of exceedance at points over some fixed time and frequency, usually for building code applications. Time‐dependent probability calculations are based on the concept of elastic rebound and assume rupture quasi‐periodicity that leads to higher probability associated with a longer elapsed time since the most recent rupture of similar magnitude from a given source. Applications are usually financial in nature and tend to be applied over shorter periods than the lifespan of a building. Operational earthquake forecasting is mostly concerned with aftershocks and can be produced in near‐real time following a mainshock or foreshock event. These short‐term forecasts are usually based on empirical time and space relations, or stress change and rate‐state models.

The possibility that aftershock faults exist complicates seismic hazard modeling in a few ways. For standard PSHA, these faults are not independent sources and are not likely to rupture randomly. Instead, an aftershock fault will only rupture after being loaded by slip on adjacent faults, and perhaps only after having such interactions being repeated. However, in time‐independent PSHA, the largest events in a spatial zone can be treated as mainshocks, regardless of their provenance (e.g., Boyd, 2012).

The impact of aftershock faults is greater for time‐dependent earthquake probability calculations. Following interaction loading, the aftershock fault would obey an Omori law decay in rupture probability versus time (e.g., Omori, 1894; Dieterich, 1994). This behavior is time‐dependent but inconsistent with elastic rebound time‐dependent models, because there is no direct monotonic tectonic loading. Thus, the idea of increasing probability with time is, in fact, reversed for aftershock faults. Aftershock faults could be appended to elastic rebound‐based time‐dependent analyses by linking them to causal ruptures and then having them follow an Omori decay model. This sort of approach is already enabled in the UCERF3‐ETAS model (Field et al., 2017). Aftershock faults are of course consistent with operational earthquake forecasts, because Omori law time dependence is their primary temporal basis.

Complex multisegment earthquake ruptures are observed globally (e.g., Eberhart‐Phillips et al., 2003; Hamling et al., 2017), and these types of ruptures are increasingly considered in seismic hazard analysis (e.g., Page, 2021) and were included in the UCERF3 effort (Field et al., 2014). Many (46%) of our identified aftershock fault subsections are included in the UCERF3 rupture set, as rupturing along with subsections of the San Andreas, Elsinore, San Jacinto, and San Gregorio faults within the San Andreas system. Simultaneous rupture of aftershock faults linked with the San Andreas system is consistent with $time−1$ Omori law time dependence, because the most likely occurrence is at near‐zero time. However, substantial delays are also observed; there does not appear to be evidence for significant slip on the Garlock fault during or after the adjacent 1857 Ft. Tejon earthquake on the San Andreas fault (McGill and Sieh, 1991; Zielke et al., 2012), despite coseismic stress transfer loading (Parsons, 2006). Similarly, the eastern Garlock fault did not have a significant response to stress loading following the 2019 M 7.1 Ridgecrest earthquake (e.g., Mancini et al., 2020), though a creep episode was observed (Ross et al., 2019). Loading of the Garlock fault only through stress interactions might explain its irregular timing of paleoearthquakes (Dawson et al., 2003).

A valid question is that how important are aftershock faults? Our preliminary analysis suggests that about 21% of all California fault subsections can be categorized as aftershock faults, though they account for about 4% of the total moment release on known faults. They do include faults that participate in ruptures with quantified 30 yr M 7.7 probability (ranging from 0.001% to 1%; Field et al., 2015), including the Anacapa‐Dume, Big Pine, Garlock, Great Valley, Hollywood, King range, Malibu Coast, Monta Vista‐Shannon, Pinto Mountain, Pt. Reyes, Santa Cruz Island, Santa Rosa Island, Santa Ynez, and Ventura‐Pitas Point faults. Figure 5 shows a distribution of aftershock fault subsection 30 yr probabilities for M 6.7 earthquakes, which can get as high as ∼5% (Field et al., 2015). According to our model concept, these probability values are only relevant if the aftershock faults rupture within multifault earthquakes controlled by faults loaded by tectonic motions.

We observe complications associated with creeping faults. For example, we identify the Coalinga section of the Great Valley fault system as an aftershock fault. This fault section ruptured the M 6.7 Coalinga earthquake in 1983. We find that San Andreas fault slip loads the Coalinga section (Fig. 3b), though this has occurred by creep rather than San Andreas fault earthquakes; Simpson et al. (1988) and Toda and Stein (2002) observed stress interactions between the creeping section of the San Andreas fault and the Coalinga section of the Great Valley system. Thus, this example of an aftershock fault is loaded monotonically by creep in a manner more akin to tectonic loading.

We find that the principal fault within the California system—the San Andreas fault—plays an important role in loading subsidiary faults. Some of these faults are also stressed by relative PAC‐SN motion and combine to release about 18% of the total moment on California faults. We identify an additional group that appears to only be loaded by interactions with the San Andreas fault and not by relative plate motions that we term “aftershock faults,” because, globally, faults that receive a sudden stress load from an adjacent earthquake have Omori law rupture rates (e.g., Parsons, 2002). These faults account for ∼4% of the total moment on California faults. We suggest that these types of faults might be handled in a separate category in time‐dependent probability calculations, because they likely rupture neither according to a random‐in‐time model nor an elastic rebound–renewal model. Our analysis is not comprehensive in that we only examined interactions between the San Andreas and surrounding faults, whereas there are many other possible interactions within the California fault system. An earthquake simulator model might be the most appropriate way to fully quantify the full suite of interaction‐loading effects.

Fault subsection data, including locations, strikes, dips, rakes, and slip rates are all available at http://wgcep.org (last accessed September 2022), as are subsection time‐dependent probability values. Stress change calculations were made using the DLC1.2.1 package written by R. W. Simpson, U.S. Geological Survey (USGS; available upon request). Finite‐element calculations were made using ANSYS commercial software. The supplemental material for this article includes tables and figures that show model evaluation efforts. Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. Government.

The authors acknowledge that there are no conflicts of interest recorded.

The authors thank Jeanne Hardebeck (U.S. Geological Survey [USGS]), two anonymous journal reviewers, and funding from U.S. Geological Survey Coastal and Marine Hazards and Resources Program.