## Abstract

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.

## Introduction

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.

## Data

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.

## Methods

### 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

*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 $\mu \u2019=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 $\Delta CS$ on each of the two nodal planes of the representative mechanism, and the stress change is considered positive if $\Delta CS>0$ on at least one plane.

To test the robustness of the results with respect to modeling assumptions, we compute $\Delta 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.

*et al.*, 1997) gives the expected rate of aftershocks $Raft$ at time

*t*after a stress step $\Delta \tau $:

*et al.*, 2020). The stress step $\Delta \tau $ is assumed to be replaceable by $\Delta 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:

### 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 $r\u22122$, and far‐field body waves decay as $r\u22121$. 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

*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

*i*and

*j*:

*i*performs better. The confidence intervals of IG are computed following Rhoades

*et al.*(2011).

## Results

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 ($<10\u22125\u2009\u2009events/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 ($\u226510\u22122\u2009\u2009events/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 ($\u226410\u22124\u2009\u2009events/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 $|\Delta CS|$ for both the stress shadow ($\Delta CS<0$) and stress increase ($\Delta 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 $|\Delta CS|$. We observe, in contrast, that the earthquake rate in the shadows generally increases, and that $Raft/Rback$ increases with $|\Delta 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/Rback\u221dr\u22122.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 $r\u22122$ is still an acceptable fit to the observations (Fig. 3e). A decay of $r\u22122$ 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).

## Discussion

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 $r\u22122$. 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 ($\u226510\u22122\u2009\u2009events/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.

## Conclusions

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 $r\u22122$, 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.

## Data and Resources

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).

## Declaration of Competing Interests

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

## Acknowledgments

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.