For decades there has been a debate about the relative effects of dynamic versus static stress triggering of aftershocks. According to the static Coulomb stress change hypothesis, aftershocks should not occur in stress shadows—regions where static Coulomb stress has been reduced. We show that static stress shadows substantially influence aftershock occurrence following three M ≥ 7 California mainshocks. Within the modeled static Coulomb stress shadows, the aftershock rate is an order of magnitude lower than in the modeled increase regions. However, the earthquake rate in the stress shadows does not decrease below the background rate, as predicted by Coulomb stress change models. Aftershocks in the stress shadows exhibit different spatial–temporal characteristics from aftershocks in the stress increase regions. The aftershock rate in the stress shadows decays as a power law with distance from the mainshock, consistent with a simple model of dynamic stress triggering. These aftershocks begin with a burst of activity during the first few days after the mainshock, also consistent with dynamic stress triggering. Our interpretation is that aftershock sequences are the combined result of static and dynamic stress triggering, with an estimated ∼34% of aftershocks due to dynamic triggering and ∼66% due to static triggering.

The triggering of aftershocks by a mainshock is frequently explained by static stress changes that permanently alter the stress field. The widely used static Coulomb stress change model (Reasenberg and Simpson, 1992) also predicts the existence of stress shadows, regions where a decrease of Coulomb stress on active faults means that no or few aftershocks should occur. Suppression of moderate‐to‐large earthquakes by a large mainshock has been observed within modeled stress shadows (Harris and Simpson, 1996, 1998). However, statistically significant rate decreases of smaller earthquakes are not found after a mainshock (Felzer and Brodsky, 2005).

The occurrence of some aftershocks inside stress shadows has been attributed to heterogeneity in background stress or fault orientations, making some faults anomalously favorable for aftershock production (Marsan, 2006). Alternatively, these aftershocks could be triggered by transient dynamic stresses from passing seismic waves. Dynamic stresses trigger distant earthquakes (Hill et al., 1993), and projection of dynamic far‐field triggering into the near‐field suggests that 15%–60% of aftershocks are dynamically triggered (van der Elst and Brodsky, 2010). Dynamic stresses also affect the propagation of earthquake rupture. However, the relative importance of dynamic versus static stress triggering of aftershocks in the near field, within about two fault lengths of the mainshock, remains controversial (Felzer and Brodsky, 2006; Richards‐Dinger et al., 2010).

To quantify the effects of stress shadows on aftershock triggering, we investigate the rate of aftershocks inside and outside the modeled stress shadows of three M 7.1–7.3 California mainshocks. We examine whether the spatial and temporal characteristics of these two sets of aftershocks are consistent with static or dynamic triggering.

Earthquakes before and after the 1992 M 7.3 Landers, 1999 M 7.1 Hector Mine, and 2019 M 7.1 Ridgecrest, California, mainshocks are obtained from the Southern California Seismic Network (see Data and Resources). We obtain Ridgecrest aftershocks through the time of this study (1125 days), along with similar‐length (1200 day) aftershock catalogs for the other mainshocks and 5 yr of background seismicity. The catalog in the study areas appears to be complete at M 2.0 (Helmstetter et al., 2006), although it is less complete immediately following each mainshock. To study the spatial density of aftershocks, we use M ≥ 2.0 events within the first three months after the mainshock, long enough to have enough aftershocks to study but not so long to become dominated by background earthquakes. The temporal evolution of the aftershocks is determined from the M ≥ 3.0 events, as the catalog should be complete down to this magnitude after the first day (Helmstetter et al., 2006), allowing study of the early part of the aftershock sequence.

Aftershock rate

To study the spatial distribution of aftershock rate, we compute the earthquake rate change, Raft/Rback, as a function of location, in which Raft is the aftershock rate, and Rback is the background earthquake rate. The region around each mainshock is divided laterally into a grid of 2 km × 2 km cells, and the number of background earthquakes and aftershocks occurring in each cell is recorded. We do not attempt to compute rate changes on this high‐resolution spatial grid, because of the difficulties of detecting rate decreases in the presence of low seismicity rates (Marsan and Nalbant, 2005). Instead, we partition the grid cells into 30 equal‐area sets according to their values of static Coulomb stress change or distance from the mainshock, and compute Raft/Rback for each set. The uncertainty in Raft/Rback is computed assuming Poisson distributed uncertainty for the number of aftershocks and background earthquakes. The dependence of Raft/Rback on static Coulomb stress change and distance from the mainshock is then compared to the predictions of the static and dynamic stress triggering models.

To study the temporal evolution of the aftershock rate, we separately plot the cumulative number of aftershocks Naft, inside and outside the modeled stress shadows, as a function of time after the mainshock. Secondary aftershocks triggering is removed from the catalog using the statistical epidemic‐type aftershock sequence (ETAS) model (Ogata, 1988), with parameters appropriate for California seismicity (Field et al., 2017). The predicted number of secondary aftershocks from the ETAS model is computed as a function of time, and this is subtracted from the number of observed aftershocks. We exclude aftershocks below the time‐dependent magnitude of completeness (Helmstetter et al., 2006) from the catalog and from models compared to the catalog.

We analyze grid cells ≥5 km distance from the mainshock fault planes and with the maximum shear stress change ≥0.01 MPa. The results of this work, therefore, apply to this near‐field region but exclude very near field aftershocks within 5 km of the mainshock rupture, which make up ∼60% of all aftershocks.

Static stress triggering model

The static Coulomb stress change ΔCS is defined as follows (Reasenberg and Simpson, 1992):
(1)ΔCS=Δτ+μΔσ,
in which Δτ is the shear stress change resolved in the slip direction of the receiver fault, Δσ is the normal stress change (tension positive), and μ is the effective coefficient of friction. We compute ΔCS from finite‐source slip distributions for each of the three mainshocks (Wald and Heaton, 1994; Ji et al., 2002; Liu et al., 2019) using dislocations in an elastic half space and the subroutines of Okada (1992). We include the slip distribution of the M 6.4 Ridgecrest foreshock as part of the Ridgecrest mainshock, and the M 6.1 Joshua Tree foreshock and M 6.3 Big Bear aftershock as part of the Landers mainshock. The shear modulus is assumed to be 30 GPa, and μ=0.4 (Stein et al., 1997). Stress changes are computed on the same 2 km × 2 km grid used to compute the earthquake rate at a depth of 6 km—the approximate median earthquake depth in the study areas.

Receiver faults are chosen based on the focal mechanisms of the premainshock earthquakes from the updated Yang et al. (2012) catalog. For each grid cell, we consider all quality A and B focal mechanisms of events occurring within the cell, or the closest 30 events if there are not that many within the cell. The median radius is ∼30 km. We define the representative focal mechanism for the grid cell as the average focal mechanism and define its uncertainty as the root mean square angular difference of each mechanism from this average. We compute ΔCS on each of the two nodal planes of the representative mechanism, and the stress change is considered positive if ΔCS>0 on at least one plane.

To test the robustness of the results with respect to modeling assumptions, we compute ΔCS for alternate mainshock slip models (Hernandez et al., 1999; Kaverina et al., 2002; Barnhart et al., 2019), for values of the effective coefficient of friction ranging from 0 to 0.8 in 0.2 increments, and at 2 km and 10 km depth. We also test different receiver fault orientations by randomly perturbing the representative mechanism for each grid cell based on its uncertainty.

The Coulomb Rate and State (CRS) model (Stein et al., 1997) gives the expected rate of aftershocks Raft at time t after a stress step Δτ:
(2)Raft=Rback[exp(ΔτAσ)1]exp[tta]+1,
in which Rback is the background rate, Aσ is a constant times the normal stress, and ta=Aσ/τ˙, in which τ˙ is the stressing rate. We assume Aσ=0.02  MPa, τ˙=5.5×104  MPa/  yr for Landers and Hector Mine (Gross and Kisslinger, 1997), and τ˙=2.7×104  MPa/  yr for Ridgecrest (Mancini et al., 2020). The stress step Δτ is assumed to be replaceable by ΔCS (Harris and Simpson, 1998). To find the expected number of aftershocks Naft from time t1=0 to time t2, we integrate over this time period:
(3)Naft=Rback0t21[exp(ΔCSAσ)1]exp[tta]+1dt,
which is
(4)Naft=Rbackta[t2ta+ΔCSAσ+ln([exp(ΔCSAσ)1]exp[t2ta]+1)],
and the rate change over the time period 0 to t2 is therefore
(5)RaftRback=Naftt2Rback=1+tat2[ΔCSAσ+ln([exp(ΔCSAσ)1]×exp[t2ta]+1)].

Dynamic stress triggering model

There is no widely accepted quantitative model for the rate of dynamically triggered aftershocks. The aftershock rate should be related to the dynamic strain, which roughly scales with the amplitude of seismic waves (Gomberg, 1996; van der Elst and Brodsky, 2010). Near‐field body waves decay with distance as r2, and far‐field body waves decay as r1. Modeling dynamic triggering is complicated by many factors, however, including the radiation pattern and directivity of the mainshock (Kilb et al., 2000), seismic velocity structure and attenuation, and the frequency dependence of triggering (Gomberg, 1996). We assume a simple model in which Raft/Rback decays as a power law with distance r from the mainshock fault, with an exponent fit to the observations.

Model evaluation

Models can be compared using the likelihood of the observed aftershocks, given each model (Rhoades et al., 2011). We test the spatial distribution of aftershocks, not the number; so the spatial density is normalized such that the total number of forecast events equals the observed number Nobs. The number of observed events in each grid cell is denoted by nobs(x), and the forecast number from model i is ni(x). The log likelihood of the model, summed over all grid cells x in the study area A, is
(6)LLi=xA[ni(x)+nobs(x)ln(ni(x))ln(nobs(x)!)].
We use the information gain (IG) to compare models i and j:
(7)IGij=(LLiLLj)/Nobs.
Positive values indicate that model i performs better. The confidence intervals of IG are computed following Rhoades et al. (2011).

Figure 1 illustrates the problem of aftershocks occurring where they preferentially should not within modeled static Coulomb stress shadows. The CRS model exhibits stress shadows with very low aftershock density (<105  events/km2), with the lowest density closest to the mainshock fault plane (Fig. 1a), as would be expected, because this is where the calculated stress decreases are the largest. In contrast, the smoothed spatial density of M ≥ 2 aftershocks is relatively high (102  events/km2) everywhere near the mainshock fault plane (Fig. 1b). Synthetic aftershock catalogs are generated by randomly selecting locations from the CRS spatial density model (Fig. 1a) and applying the same adaptive smoothing method used on the real data. The smoothed synthetic catalogs exhibit stress shadows with very low seismicity rates (104  events/km2) near the mainshock (Fig. 2), so the smoothing cannot explain the lack of stress shadows in the real catalog.

We plot Raft/Rback as a function of |ΔCS| for both the stress shadow (ΔCS<0) and stress increase (ΔCS>0) regions (Fig. 3a–c). The stress shadows appear to have a substantial effect on the earthquake rate change. The observed Raft/Rback in the modeled stress shadows is approximately an order of magnitude lower than in the stress increase regions, for comparable absolute Coulomb stress changes greater than about 0.01 MPa. The CRS model fits reasonably well for the stress increase regions but not for the stress shadows. The CRS model predicts earthquake rate to decrease to below background in the shadows, with Raft/Rback rapidly decreasing with increasing |ΔCS|. We observe, in contrast, that the earthquake rate in the shadows generally increases, and that Raft/Rback increases with |ΔCS|.

We next investigate whether dynamic triggering can explain the rate of aftershocks in the stress shadows (Fig. 3d–f). We fit a power law to the observed Raft/Rback versus distance r, out to ∼80 km at which distance Raft/Rback becomes indistinguishable from 1. We find that Raft/Rbackr2.0 is a good fit to the rate changes in the stress shadows. For the Landers and Ridgecrest sequences, the best‐fitting exponent is −2.0 and −1.9, respectively. For the Hector Mine sequence, the best‐fit decay is somewhat steeper with an exponent of −2.3, whereas a decay with r2 is still an acceptable fit to the observations (Fig. 3e). A decay of r2 is consistent with the decay of near‐field body waves. Alternatively, it may be consistent with far‐field body waves experiencing attenuation. Attenuation is an exponential function of travel time (or equivalently, distance), rather than a power law, but the difference may not be discernable given the uncertainty and scatter of the observations. Other studies have found a similar decay for full catalogs (Felzer and Brodsky, 2006).

We next examine the temporal evolution of aftershock rate. In the Coulomb stress increase regions, the aftershock rate is fit well by the CRS model, which approximates a power law decay with time (Fig. 4a–c). The aftershocks in the stress shadows, in contrast, begin with an initial burst during the first 6–13 days, followed by additional aftershocks occurring at a much slower, approximately, constant rate (Fig. 4d–f). (The Landers sequence exhibits some additional activity about 150 days after the mainshock and then returns to the same slower rate.) An initial burst of aftershocks is consistent with dynamic stress triggering, which is the most important early in the aftershock sequence (Marsan and Nalbant, 2005). Dynamic stress changes can directly trigger earthquakes as the seismic waves pass and can cause delayed triggering through perturbations to properties such as pore pressure or crack growth (Brodsky and van der Elst, 2014). Most of the aftershocks within the first 90 days occur during the initial burst, so it is the spatial characteristics of this burst that are represented in Figure 3. The subsequent slow, apparently constant rate of aftershocks may reflect the beginning of a much longer timescale process such as viscoelastic relaxation (Pollitz and Cattania, 2017).

Our results suggest that aftershock triggering is the combined result of two physical processes: Coulomb static stress triggering that is modeled reasonably well by CRS and dynamic stress triggering that produces a burst of early aftershocks that decay with distance from the mainshock as approximately r2. In Coulomb static stress shadows, in which static stress changes generate few or no aftershocks, dynamic triggering is the dominant mechanism. In the static stress increase regions, both the mechanisms are active.

We estimate the fraction of aftershocks that are dynamically triggered in this hybrid static–dynamic stress triggering model. We assume that all aftershocks in the stress shadows during the first three months are dynamically triggered, and that a similar rate of dynamically triggered aftershocks exists within the positive Coulomb stress regions. This may be an underestimate, as the dynamic stress change may be higher in lobes of positive static stress (Kilb et al., 2000), meaning that the rate of dynamically triggered aftershocks may be higher in positive static stress change lobes than in negative stress change lobes. The computed fraction of dynamically triggered aftershocks should therefore be considered a lower bound. Aftershocks can trigger secondary aftershocks, usually close by, which would affect the apparent rates of static‐ and dynamic‐triggered aftershocks proportionally. To test the robustness of our estimate with respect to modeling assumptions, we repeat this calculation for alternate mainshock slip models, multiple receiver fault orientations, and different values of the coefficient of friction. We find a median value of 34% of aftershocks dynamically triggered, with 68% of the 600 realizations falling between 26% and 48% (Fig. 5a).

A hybrid model of aftershock density that weights the CRS model 66% and the simple dynamic triggering model 34% appears much more similar to the observed aftershock density than the CRS‐only model (Fig. 1c). Similar to the observed aftershocks, the hybrid model exhibits relatively high aftershock density (102  events/km2) everywhere near the mainshock fault plane. The hybrid model is still missing some features of the observed aftershocks, for example, the directivity that causes higher aftershock rate to the north of the Landers rupture (Kilb et al., 2000). Other possible triggering mechanisms such as postseismic relaxation, slow afterslip, and pore pressure change are also not included in the hybrid model.

The hybrid model performs quantitatively better in forecasting the aftershock locations than the CRS model alone (Fig. 5b). The hybrid model generally performs better in forecasting than the dynamic model alone for the Landers and Hector Mine sequences, whereas the dynamic model performs better for the Ridgecrest sequence (Fig. 5c). The dynamic‐only model is similar to empirical aftershock models based on distance from the mainshock (Felzer and Brodsky, 2006). Our results imply that the hybrid model is often, but not always, an improvement over such models.

The hybrid static–dynamic stress triggering model can reconcile some apparently conflicting previous results. Dynamic stresses may cause numerous small aftershocks to occur inside static stress shadows, so that the theoretical decrease in the rate of earthquakes due to static stress decrease is not detectable (Felzer and Brodsky, 2005). On the other hand, the order of magnitude greater rate of aftershocks in the stress increase lobes could explain why static stress shadows are visible in other studies (Harris and Simpson, 1996, 1998). The order of magnitude difference in rate also suggests that static stress shadows may be useful in forecasting the spatial distribution of aftershocks.

The spatial–temporal characteristics of aftershocks in the stress shadows are better explained by the hybrid static–dynamic stress triggering model than by CRS triggering on faults of anomalous orientation. If CRS were the main triggering process in both the stress increase and stress shadow regions, the aftershocks in the two regions should have similar spatial and temporal decay. The main difference would be a lower aftershock rate in the stress shadows, scaling with the lower rate of background earthquakes on anomalously oriented faults. Instead, the aftershocks in the two regions have different forms of both the decay with distance from the mainshock and the time evolution. These differences imply different dominant triggering processes.

Aftershock triggering in the very near field is a difficult problem that we do not attempt to address in this study, so we have excluded the region within 5 km of the mainshock from our analysis. Static stress changes in the very near field are difficult to calculate, because mainshock slip models usually have resolution of at the best a few kilometers. In addition, static stress changes in the very near field can be highly dependent on small‐scale fault geometry features. Many static stress triggering studies compute Coulomb stress changes on optimally oriented fault planes (Stein et al., 1997), which eliminates the modeled stress shadow very close to the fault. These modeling issues obscure the relative importance of dynamic versus static stress triggering in the very near field.

We find that aftershocks in the stress shadows of three M > 7 California earthquakes exhibit different spatial–temporal characteristics than aftershocks in the stress increase regions. Aftershocks in the stress increase regions are well explained by the CRS model. The aftershocks in the stress shadows occur in a burst during the first few days, and decay as a power law with distance from the mainshock as approximately r2, consistent with a simple model of dynamic stress triggering. We propose that aftershock triggering is the combined result of Coulomb static stress triggering and dynamic stress triggering. Static stress triggering dominates in the stress increase regions, whereas dynamic triggering dominates in the stress shadows. We estimate that ∼34% of aftershocks during the first three months are due to dynamic triggering and ∼66% due to static triggering.

Dynamic triggering appears to be an important contributor to aftershock triggering in the near field. However, most studies of near‐field aftershock triggering focus on static stress changes alone. A better understanding of near‐field dynamic triggering would support development of a quantitative spatial–temporal near‐field dynamic triggering model with the same level of sophistication as the CRS model for static stress triggering. The development of hybrid static–dynamic triggering models should greatly improve our ability to forecast the locations of near‐field aftershocks.

Earthquake data were obtained from the Southern California Seismic Network (SCSN), doi: 10.7914/SN/CI, operated by the Caltech Seismological Laboratory and the U.S. Geological Survey (USGS). The earthquake catalog was accessed through the Southern California Earthquake Data Center (SCEDC), doi: 10.7909/C3WD3xH1, https://service.scedc.caltech.edu/eq-catalogs/date_mag_loc.php (last accessed August 2022), and the focal mechanisms (Yang et al., 2012) were accessed at https://scedc.caltech.edu/data/alt-2011-yang-hauksson-shearer.html (last accessed February 2022). All other data came from published sources listed in the references. The Okada (1992) code to model dislocations in an elastic half space is publicly available at https://www.bosai.go.jp/e/dc3d.html (last accessed October 2021).

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

The authors are grateful to Joan Gomberg, Nicholas van der Elst, Debi Kilb, and an anonymous reviewer for their helpful comments on the article.

Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. Government.

This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.