Dynamic Earthquake Triggering in Southern California in High Resolution: Intensity, Time Decay, and Regional Variability

Earthquake triggering by seismic waves has been recognized as a phenomenon for nearly 30 years. However, our ability to study dynamic triggering has been limited by our ability to capture the triggering stresses accurately and record the resultant earthquakes. Here we use full waveforms from a dense seismic network and a modern, high‐resolution seismic catalog to measure triggering in Southern California from 2008 to 2017 based on interevent time ratios. We find that the fractional seismicity rate change, which we term triggering intensity or triggerability, as a function of peak strain change for the period of ∼20 s due to distant earthquakes is monotonically increasing and compatible with earlier measurements made with a disjoint data set from 1984 to 2008. A triggering strain of 1 microstrain is equivalent to the local productivity generated by an M1.8 earthquakes. This result implies that a prediction of seismicity rate changes can be made based on recorded ground shaking using the same formalism as currently used for aftershock prediction. For a teleseismic event, this small level of triggering occurs throughout the region and thus aggregates to a regional effect. We find that the triggering rate decays after the triggerer follows an Omori‐Utsu law, but at a much slower rate than a typical aftershock sequence. The slow decay rate suggests that an ancillary process such as creep or fluid flow must be part of dynamic triggering. The prevalence of triggering in areas of creep or fluid involvement reinforces this inference. A triggering cascade of secondary earthquakes is insufficient to explain the data.

The timing of the triggering is particularly important to quantify. Dynamically triggered earthquakes often occur subsequent to the passage of the perturbing stress in the seismic waves. This delay is hard to understand in a simple Coulomb failure model and is a fundamental feature of dynamic triggering that makes it a physically mystifying process (Belardinelli et al., 2003;Gomberg, 2001;Parsons, 2005). Capturing and quantifying the timing of the triggered events is one of the major observational gaps in dynamic triggering (Brodsky & van der Elst, 2014). However, datasets that can approach this problem have been few (Brodsky, 2006;Shelly et al., 2011).
Our ability to detect earthquakes is primarily limited by the magnitude of completeness of extent catalogs. Establishing statistically significant triggering rates requires measuring a large population of earthquakes. Since many more earthquakes are small than large, decreasing the magnitude of completeness can dramatically increase the visibility of earthquake rate changes. Recent progress in catalog generation by template matching technique has now resulted in catalogs that are much more amenable to dynamic triggering studies. The 10-fold increase in detected seismicity in the recently released Quake Template Matching (QTM) earthquake catalog of Southern California allows a similar increase in the detectability of rate changes (Ross et al., 2019). The more robust measurement of rate changes allows a more refined view of dynamic triggering to capture issues such as timing.
Measuring the triggering stresses has also previously been limited by the available seismograms. Studies opted either to confine their attention to areas with high-quality instrumentation (Gomberg et al., 2004;Miyazawa, 2019;Miyazawa & Brodsky, 2008;Velasco et al., 2008) or to extrapolate attenuation relationships to infer ground motion far away from the instruments (van der Elst & Brodsky, 2010). The attenuation relationship approach has been successful in allowing broad-scale studies but requires a consistent magnitude determination procedure, which is problematic for composite catalogs. This issue is brought into particular focus by the changes in the Advanced National Seismic System (ANSS) catalog since the publication of van der Elst and Brodsky (2010). Although, the original method still works on the original data as downloaded in 2009, the reported magnitudes for the global events in the study period of 1984-2008 have changed significantly enough since then that the method no longer works on the ANSS catalog as downloaded today. The source of the discrepancy is likely the changing magnitude types, which would require a different set of attenuation relationships to mimic the triggering strains inferred from the old magnitude decisions. A much robust approach is to use the waveforms directly to measure the triggering strains. Such an approach is now possible due to the increase in publicly available seismic records and, equally importantly, an increase in the computational capability required to process them. Direct measurement of strains is also important for rare, large magnitude events where individual variations in radiation patterns can make local amplitudes deviate strongly from the average attenuation curve inferred based on magnitude and distance.
This study uses the waveforms of Southern California and the QTM catalog to measure dynamic triggering in 2008-2017. After reviewing the mathematical framework for measuring triggering, we first establish that seismicity rate changes are related to peak dynamic strain as measured as previously proposed on an entirely disjoint data set. We then proceed to use this expanded data to measure the timing and spatial variation of triggering. Fully utilizing this data requires the development of a new statistical model for the rate changes, which is shown to match the observational data well. We find that the dynamically triggered seismicity decays more slowly with time than conventional aftershock sequences. The observation suggests that dynamic triggering involves a distinct process such as triggered creep or fluid movement.

Quantifying the Triggering Intensity
We intend to study the fractional rate changes suddenly from an initial, Poissonian background rate λ 1 to a new one λ 2 at the time of the triggering wave, then the triggering intensity or triggerability, n, can be defined as rates can be either constants or functions of time that are being estimated directly before and after the triggering waves.
An estimate of the rate change is model dependent as a Poisson process rate is not a directly observable quantity. We choose here to connect these model-dependent quantities to the observed interevent times using the interevent time ratio R as where t 1 and t 2 are times to and from the arrival of triggering waves ( Figure 1). This time ratio metric is useful as a robust way to study the statistical behavior of a large population of earthquakes and can be rigorously associated with rate changes under certain simplifying assumptions. For instance, for a homogeneous Poisson process with a step-change in seismicity rate λ due to the triggerer, i.e., the triggering wave, van der Elst and Brodsky (2010) showed that the expectation of R is By numerically solving Equation 3, we can obtain the triggering intensity n. Note that this approach presupposes a good estimate of 〈R〉, which means that an adequately sized population needs to be measured.
An alternative model for the rate change is that the rate after the triggering stress is a function of time. For instance, it might be reasonable to assume that the rate change follows a systematic decay behavior such as an Omori-Utsu law. This situation would occur if the triggering waves immediately trigger some local events during their passage and the subsequent events are aftershocks of the initial, waveform triggered earthquakes (Brodsky, 2006). The mathematical statement would be where the background rate λ is the assumed to be equal to λ 1 and K, c, and p are the model parameters.
MIYAZAWA ET AL. In Appendix A, we develop the relationship between the distributions of t 1 and t 2 under this model and find that the ratio of the cumulative number of first events from and to the triggerers is obtained as for p < 1 and for p < 1 and small t in particular. The cumulative number ratio is inversely proportional to the time to the power of p for small t. Note that this derivation is based on comparing the independent distributions of t 1 and t 2 without consideration of their relative values for a particular event. It is therefore distinct from the thinking behind Equation 2, where the ratio of t 1 and t 2 is considered for a single triggerer.
The aftershock productivity K generally increases with the magnitude of the mainshock. Another well-known seismicity model includes the epidemic-type aftershock sequence (ETAS) model (e.g., Ogata, 1988), where the seismicity rate is represented as where i M is the magnitude of the earthquake at time t i , th M is the lower threshold of magnitude, and C and α are the constants. Assuming that, the seismicity rates of both the Omori-Utsu and ETAS models are the same at the time one day after the mainshock, i.e.,    1 i t t c , there are few contributions from other large earthquakes before the mainshock, and α = 1, we have the magnitude corresponding to K as   th log .
We call this magnitude, E M , as an equivalent mainshock magnitude to the local K. The mainshock magnitude is equivalent in the sense that a local earthquake of this magnitude would trigger as many earthquakes in one bin at the time one day after the mainshock as the observed dynamic triggering.

Seismicity and Seismic Data
We implement this method by dividing the region of Southern California into spatial bins of 0.1° × 0.1° as in the previous studies (Brodsky & van der Elst, 2014;van der Elst & Brodsky, 2010). Here we use local earthquakes with M ≥ 0.5, considering the completeness magnitude, at the depths from 0 to 15 km in the QTM earthquake catalog (2008-2017) with a minimum detection threshold of 9.5 × MAD (Ross et al., 2019), whereas the previous studies used local earthquakes of M ≥ 1.5 from the Southern California Seismic Network (SCSN) catalog. The individual R for each triggerer is calculated in each bin when the completeness magnitude M c ≤ 0.5 ( Figure 2). The magnitude-time diagram of seismicity for the first event to/from the triggerer clearly shows the increase of seismicity immediately after the triggerer within one day for large peak ground velocities (  The peak amplitude of seismic waves from a distant earthquake is directly obtained from the waveforms recorded at SCSN broadband seismic stations. We first select global earthquakes with M ≥ 5.0 from ANSS Comprehensive Earthquake Catalog during the term from 2008 to 2017, which corresponds to the period that the local earthquake catalog is available (The previous study used earthquakes with M ≥ 5.0 from 1981 to April 2008). Since we are interested in the dynamic triggering, the distance between a triggering earthquake and a local spatial bin is set to be greater than 50 km. We also use the largest earthquake in terms of magnitude in any 2-hour time window to avoid redundantly picking a triggering peak in the continuous waveform. The number of triggering events thus selected is 11,346 out of 18,957 in total. For each triggering event, we extract continuous records with the time-window from 500 s before the theoretical P-wave arrival and to 500 s after the surface wave arrival with an apparent velocity of 3.4 km/s. The peak ground velocity (PGV) is the maximum amplitude of the vector velocity from the three components of a triggering event, band-pass filtered between 0.02 to 0.10 Hz, recorded at each seismic station with 1 Hz sampling rate. At these long periods, site effects are relatively small. A representative PGV at each bin from each event is given by the median value of PGVs at three stations closest to the bin center, after removing outliers, that is, values larger than 9.0 × MAD, which is required for excluding artificial noise. The corresponding peak dynamic strain is approximately given by dividing the PGV by a seismic wave velocity of 3.5 km/s. Note that we implicitly assume the triggerer is the largest amplitude of surface waves with the predominant period of about 20 s. The peak strain smaller than 2 × 10 −9 is not used because it is too small to trigger and is comparable to the background noise level. This means that distant earthquakes with magnitudes around 5.0 are too small to cause PGVs that exceed this strain level based on the empirical prediction, however, the predicted PGVs may underestimate the measured ones because of neglecting radiation pattern effects.
MIYAZAWA ET AL.  Therefore, those waveforms are systematically used in the initial step mentioned above. Thus obtained peak strain changes as a function of arrival time from 2008 to 2017 are shown in Figure S1, where the number of earthquakes is 4,542. The time of the observed PGV is used to measure t 1 and t 2 in each bin, while previous studies used the origin time of triggerer. The local earthquakes in each bin in the QTM catalog are about tenfold more densely populated than in the conventional SCSN catalog and the correction for travel time of triggerer may not be negligible. The R-value is calculated when both t 1 and t 2 are smaller than two years.
There may be small local earthquakes missing from the QTM catalog during and after the passage of large-amplitude waves. Figure S2 shows the frequency-magnitude distribution in the regions of M c ≤ 0.5 within one hour after surface wave arrivals from the 2009 M w 6.9 Gulf of California earthquake, the 2010 M w 7.2 El Mayor-Cucapah earthquake, and the 2017 M w 8.2 Tehuantepec earthquake, where the peak strain changes exceed 1 × 10 −6 in Southern California ( Figure S1). For the 2009 and 2017 earthquakes, where the distance is from a few hundred km to a few thousand km, it likely holds M c ≤ 0.5. On the other hand, after the 2010 earthquake that occurred within a few hundred km from the study region, M c becomes larger than 0.5 due to the numerous aftershocks. This means that the n-value may be underestimated immediately after the 2010 earthquake. Moreover, M c two to four hours after mainshock seems to return to the background level smaller than 0.5, therefore t 2 larger than two hours would be robust even in this specific case. Figure 4 shows the triggering intensity n as a function of peak strain changes. The vertical error bar indicates a 90% confidence interval of R values examined by 1,000 bootstrap replications. The n-value for low strain change shows a small value with small errors, indicating no significant seismicity rate change associated with small strain perturbations. The error bars generally become larger for the greater strain changes due to the reduced number of R-values for each bin. Most of the n-values are positive, indicating a rate increase. A significant increase of n-value can be seen when the peak strain change is larger than 5 × 10 −7 , which is about one order of magnitude larger than the strain change associated with earth tide, where 82% and 89% bins are positive for the three components and the vertical component, respectively. The local maximum and local minimum appeared around 1 × 10 −6 and 2 × 10 −6 , respectively, for the three-component waveforms, which may be because of the sparsity of observations for these peak strain changes (see Figure S1).

Triggering Intensity in Southern California
Overall, it is clear that the triggering intensity n increases with the triggering strain, and the relationship is consistent with the exponential model obtained by the previous study (Brodsky & van der Elst, 2014). The relationship between triggering intensity n and peak microstrain ε can be fit by The log-log plot shows that the present result seems to have a slightly steeper curve than the previous one. Similar characteristics are observed when the peak strain change is obtained from the vertical component of waveforms. The depth dependence of n-values is also examined for every 5 km from 0 to 20 km ( Figure S3), where a similar relationship between the n-value and peak strain change is seen at every depth. The largest amplitude waves appear to trigger more local earthquakes at 0-5 km depth than larger depths. There might MIYAZAWA ET AL.
10.1029/2020AV000309 6 of 15  therefore be a slight increase in triggerability at shallow depths. More importantly, the repetition of the overall triggering pattern demonstrates the robustness of the fundamental result of increased triggered rate with increased peak strain change.

Delay Triggering Timing
Once the seismicity has increased, it should eventually return to the background rate; otherwise, the seismicity rate would grow without bound with repeated perturbations. We now examine how late the dynamic triggering can occur and how the triggered seismicity decreases with time by comparing the distribution of observed times t 1 to that of t 2 with the model of Equation 5 which is based on the assumption of a power-law decay of the dynamic triggering. Figure 5a shows the cumulative number of events with t 1 and t 2 for greater than a given value. For the purpose of an initial analysis, we limit the data to strong triggering events, that is, those with local strain changes more than 1 × 10 −6 . Since the time of events, t 1 and t 2 , are sorted in ascending order for statistical analysis, the association of t 1 and t 2 for a given triggerer is not preserved. The cumulative number of events at small times t 2 is generally larger than that for t 1 because the triggerer advances the occurrence of earthquakes, that is, n > 1, as is shown in Figure 4.
From these distributions, we can compute the ratio of the distribution, r o (t). In Figure 5b, we compute a value of r o for every observed value of t 2 with the corresponding value of the denominator from a linear interpolation of the t 1 distribution in Figure 5a. The ratio rapidly increases to reach the peak value around t = 0.01 day, which may indicate fewer observations or early triggering before the arrival of peak strain (Figure 3c). Then the ratio decreases gradually with time as expected. Note that there is a continuous function from the peak to the long delays of more than 10 days, suggesting that a single process accounts for this decay.
The model has three parameters: background seismicity λ, aftershock activity K, and the decay exponent p. One of the key features of the new catalog is the ability to estimate the background seismicity independently because of the lower completeness threshold (Ross et al., 2019). Therefore, we estimate the background seismicity for 30 days prior to the triggerer and then fit the data of Figure 5b that shows a rate decay with our model of Equation 5 and only two free parameters: K and p. To objectively limit the data to the decaying time and eliminate the early period that is likely incomplete due to observational limitations, we limit the data to that at least two hours after the triggerer or the time when the cumulative number for t 2 is at least 20, whichever comes later (black in Figure 5b).
The fitted model curve (dashed line in Figure 5b) shows that our model based on the Omori-Utsu law for the triggered earthquakes, successfully explains the observation. The resulting values of p and K values are 0.4077 ± 0.0122 and 0.08221 ± 0.00297, respectively, with a 90% confidence interval. The small value of the fit parameter p is significantly less than the usual p-value of about one for aftershock sequences (Utsu et al., 1995). Note that the model of Equation 5 assumed a small value c, which is a parameter that represents the time before the seismicity rate begins to decay and is typically much less than 0.1 days (e.g., Enescu et al., 2009). If this assumption of a low c-value is erroneous, it could potentially artificially depress p. We therefore relax the low c-value assumption and fit the data allowing c to vary. We find that p = 0.4151 and c = 2 × 10 −7 . The alternative fit indicates the approximation neglecting c satisfies the best-fit to the data. We also note that if the value of c is the large enough to influence the apparent fit, it would be much larger than MIYAZAWA ET AL.  for ordinary aftershocks and thus also indicate that the decay is slow. Thus, regardless of the parameterization, the physical conclusion is that the decay of dynamic triggering is slower than ordinary aftershocks.
This slow decay is a crucially important characteristic of remotely or dynamically triggered seismicity and constrains the mechanism of delayed triggering. The observable triggered seismicity returns to the background seismicity within a few days (Figure 3), although physical processes that prolong the process beyond the current observational level are still possible.
We now turn our attention to how the seismicity decay might vary for events separated by the amplitude of the local strain at the triggering site. Separating by strain avoids mixing different triggering strengths and can potentially illuminate the dependence of triggering productivity K on the imposed dynamic strain. Therefore, we separated the data for strain bins of 2.5 × 10 −7 to 5.0 × 10 −7 , 5.0 × 10 −7 to 1.0 × 10 −6 , 1.0 × 10 −6 to 2.0 × 10 −6 , 2.0 × 10 −6 to 4.0 × 10 −6 , and 4.0 × 10 −6 to 8.0 × 10 −6 and attempted to repeat the procedure leading to Figure 5b for each bin to obtain the model parameters K and p as a function of triggering strain. The corresponding distributions of t 1 , t 2, and r o are shown in Figures S4 and S5. The results show that for very small strains, the weak triggering signal posed some problems. Perhaps unsurprisingly, for strains less than 2.0 × 10 −6 , r o is close to 1 and the distributions of t 1 and t 2 are too similar to measure a refined feature of the data like the triggering productivity K. The strain bin of 1.0 × 10 −6 to 2.0 × 10 −6 has a more interpretable, but still marginal result as the cumulative number of t 1 shows non-Poissonian high activity ( Figure S4), which causes no increases in the ratio. Therefore, for this bin, we use an alternative method and obtain the parameters directly from the t 2 -curve in Figure S4 using Equation A7. In this alternative approach, there are four free parameters as opposed to the two free parameters of the preferred approach. Thus estimated values for 1.0 × 10 −6 to 2.0 × 10 −6 may be poorly constrained and not be directly compared with the other strain bins. For 2.0 × 10 −6 to 4.0 × 10 −6 and 4.0 × 10 −6 to 8.0 × 10 −6 , we obtain stable results because the data of wide time range can be used to fit the curve, where the slope for small t becomes steeper for larger strain.
The resulting values of K and p in the resolvable strain bins are in Figure 6. The productivity K-value becomes larger for higher peak strain change, which is consistent with the results seen for triggering intensity (Figure 4). The relationship between triggering productivity K and microstrain ε can be fit as a power law of the form (open circle with error bars) for a given range of peak strain change (horizontal bar), and plotted in the middle of the strain range used. The vertical error shows a 90% confidence interval. The gray horizontal bar shows the value for the peak strain change is larger than 1.0 × 10 −6 and smaller than the maximum peak strain change. Solid and dashed lines in (a) are regression fits, Equations 10 and 11, assuming a power law and a linear function, respectively. The equivalent magnitude corresponding to K is shown in the right vertical label. It is useful to interpret the level of dynamic triggering in terms of an equivalent mainshock magnitude of the triggering M E . The seismicity rate of dynamically triggered earthquakes from a given level of strain is equal to the seismicity rate of aftershocks that would be triggered one day after the mainshock by the equivalent mainshock magnitude as defined by Equation 8. We implement Equation 8 by estimating the ETAS parameters in Southern California in Equation 7 using a maximum likelihood ETAS inversion (Brodsky & Lajoie, 2013;Ogata, 1988). Since the QTM catalog is too large to solve the inverse problem all at once, we use every 10,000 events and get an average value for the aftershock productivity C of 0.0027 and that for the decay factor p of 1.06 with a standard deviation of 0.23. The equivalent mainshock magnitude for K is shown in Figure 6a. The triggering strain of one microstrain is equivalent to an earthquake with an equivalent magnitude of 1.8.
We now focus on the fit of the decay exponent p. The relationship between p-value and peak strain change is not clear. Though the slope seems to become steep for peak strain change larger than 2.0 × 10 −6 ( Figures S5  and S6), the p-value may not be constrained because we do not use the ratio for very small t. In all cases, the derived value of p is lower than 0.5, which is much below the conventional aftershock value of about 1.0. These results can be contrasted with p-values for remotely triggered seismicity obtained elsewhere by a conventional approach that did not show this significant depression. Individual remotely triggered sequences in the geothermal regions, the Geysers, CA, and Oita, Japan, had p-values of 1.03 and 1.81, respectively (Brodsky, 2006;Miyazawa, 2016). At the Coso geothermal region, we again fit p = 1.81 for triggered seismicity from the 1992 M w 7.3 Landers earthquake with respect to a peak strain change of ∼1 × 10 −5 ( Figure S7). A smaller, but still significant discrepancy exists with the results of Pankow et al. (2004), who reported p = 0.53 to 0.75 for remotely triggered seismicity at Utah within 25 days by the 2002 M w 7.9 Denali earthquake with peak strain changes of about 0.4 × 10 −5 to 1.0 × 10 −5 , which is only slightly smaller than the mean p-value of the study region of about 0.8. Other studies obtain small p-values of about 0.7 for triggered seismicity using a conventional approach in California (e.g., Aiken & Peng, 2014;Enescu et al., 2009) and western Japan (Opris et al., 2018) and interpret the depression of p-value as a result of stacked swarm-like earthquake sequences. However, even the stacked sequences do not have a p-value as low as observed here. The individual study p-values are not directly comparable to the present measurement as the model of Equation 5 is designed to capture the seismicity rate decay directly from the remote triggering as recorded by the cumulative number ratios, whereas the individual studies are measuring the decay of the seismicity in a single location. We will return to this discrepancy in the discussion section below. That said, at least one individual study did document a gradual decay. Meng and Peng (2014) found a gradual decay exponent of 0.65 for a dynamically triggered sequence in the Salton Sea Geothermal Field, but did not provide a reference value of the region.
The triggerability obtained using Equation 3 assumes a homogeneous Poisson process with a step-change in the seismicity rate due to the triggerer. Given the results above, we may choose to modify this equation for the non-stationary (Omori-Utsu decaying) Poisson process. However, we do not adopt this model because van der Elst and Brodsky (2010) showed the difference in 〈R〉 is small and much more sensitive to the seismicity rate change than the choice of decay model (See Figure 3 of van der Elst & Brodsky, 2010).
It is tempting to think that the depressed p-value is the result of mixing non-triggered with triggered regions. However, the model of Appendix A includes these effects through the explicit recognition of both background and triggered events in determining the observed number of earthquakes. Figure 7a shows the spatial distribution of triggering intensity n in Southern California. The map is smoothed by a two-dimensional Gaussian kernel with a standard deviation of 0.3°, which is the same approach as the previous study (van der Elst & Brodsky, 2010). This smoothing corresponds to an estimation of 〈R〉 statistically securing enough number of observations. For some isolated bins, which have significant values, sufficient smoothing may not be applied to compare with other values. The largest 4% of triggering values (n > 0.03) are observed at three regions; (i) San Andreas fault zone (118.0°W, 34.5°N), (ii) the vicinity of the Beta offshore platform producing oil and gas, and (iii) Coso volcanic field. The region along the San Jacinto fault zone also has a relatively large value. These high-value distributions are also seen in the triggering intensity map by van der Elst and Brodsky (2010). The high triggering near the Beta offshore platform and in the Coso volcanic field is consistent with the fact that triggered seismicity is often observed in volcanic and geothermal areas or areas of anthropogenic subsurface activities, where the hydrologic systems play an important role in triggering (e.g., Hill & Prejean, 2015;van der Elst et al., 2013). The vicinity of the Beta offshore platform may be stimulated due to the gas and oil production so that the earthquake can be sensitive to passing seismic waves. These regions with high n-values are located within several tens of km from the production platform, but such a spatial discrepancy can be explained by the stress transfer through poroelastic effects as the remotely induced seismicity is observed in other injection fields (e.g., Goebel & Brodsky, 2018;Goebel et al., 2019;Segall & Lu, 2015).

Spatial Changes in Triggerability
The n-values at the Coso volcanic field are also large. Previous research showed dynamic earthquake triggering in Coso (e.g., Aiken & Peng, 2014;Peng et al., 2010;Prejean et al., 2004;Roquemore & Simila, 1994;Zhang et al., 2017), but the individual case studies do not capture the trend as a function of seismic wave amplitude accessible from this study. Zhang et al. (2017) used a finer spatial grid to observe reduced triggered seismicity within the Coso geothermal field than in the surrounding areas. This absence of triggered seismicity may result in reducing the total triggering intensity in this region because the present study uses a larger spatial bin, and the smoothing that includes surrounding values can eventually illuminate high values.
The geothermal regions do not always correspond to large n-values. For example, the southern Salton Sea is a seismically active geothermal region; however, the n-value map does not show consistently high values. Since the number of triggered seismicity cannot be quantified in terms of n-value that is a normalized seismicity rate change, we obtain a value of n-value (Figure 7a) times the number of earthquakes N (Figure 2b) in each spatial bin, which is approximately equivalent to the number of triggered earthquakes (Figure 7b). The total number of triggered events is high at the southern Salton Sea.
There are some areas of high triggered seismicity along the faults in Figure 7b, while the value depends on the population of seismicity in Figure 2b. In these regions, many numbers of triggered seismicity were observed. The significant one is observed in the Anza section along the San Jacinto fault (116.5°W, 33.5°N), though this may be due to the extremely high background seismicity and aftershocks associated with the 2010 M w 5.4 and 2016 M w 5.2 Borrego Springs earthquakes, together with slightly high triggering intensities. The creep at a deep extension of the San Jacinto fault in the crust has been suggested (Wdowinski, 2009), and the delayed remote triggering of earthquakes from the 2014 M w 7.3 Papanoa earthquake is interpreted that the immediate triggering of creep might cause stress build-up at the locked fault to trigger the usual earthquakes (Li et al., 2019). Other high values in the region south to the Salton Sea (116.0°W, 32.8°N) may be associated with the aftershocks of the 2010 M w 7.2 El Mayor-Cucapah earthquake.
The correlation coefficient between the triggering intensity n (Figure 7a) and the total earthquake number N (Figure 2b) is estimated to be only 0.08691, whereas the p-value of 0.1852 for the statistical test is not significant, therefore we retain the null hypothesis of no relationship. This lack of correlation seems to be inconsistent with the result by van der Elst and Brodsky (2010). As described above, both maps look similar and the previous study qualitatively pointed out the possible positive correlation based on visual inspection. Here, we instead apply a quantitative approach and arrive at a different result. One possible explanation is that the correlations between triggering intensity and background activity changed between 2008 and 2017. We view this scenario as unlikely in light of the stability of the triggering relationship as a function of dynamic strain. Instead, we infer that the present result indicates that the seismicity can be enhanced by dynamic strain changes at any place regardless of the background seismicity. The current observation directly contradicts the prediction of Dieterich (1994) that triggerability should be correlated with background rate as both variables are controlled by the regional stressing rate.

Discussion
The major results of this study are: (1) larger peak dynamic strains result in systematically larger triggering rates with fractional rate changes of 0.385ε 0.944 for microstrain ε as measured directly from the seismic wavefield, (2) the triggered earthquakes decay with an exponent of 0.41, (3) the vicinity of the Beta offshore platform, San Andreas Fault Zone, and the Coso volcanic region are the most conspicuous regions of high triggerability. These areas include regions of partial creep and/or human-induced seismicity.
The first observation is a refinement of an earlier measurement of rate changes of van der Elst and Brodsky (2010) on an entirely disjoint data set utilizing a more robust and reliable method of measuring the perturbing strain. The framing in terms of the productivity constant K shows that indeed triggering increases with perturbing strain. The result now positions us to use the 0.0494ε 1.17 relationship to anticipate future seismicity based on the current ground shaking. This method could be implemented for operational forecasting analogous to current aftershock foreshocks. A straightforward implementation stems from the extension of the observations to a determination of the local mainshock magnitude equivalence in Figure 6a. At first blush, the small values of equivalent mainshock M E appear to indicate that the dynamic triggering is insignificant. However, teleseismic and regional dynamic triggering is distinct in that an entire region can experience the triggering at once, thus dramatically increasing the total number of triggered events. For instance, if all of the study regions experienced a strain of 2 × 10 −6 , the resulting seismicity that day is equivalent to a magnitude 2.2 in each of the 234 bins. According to the QTM earthquake catalog, Southern California, on average, has 0.86 magnitude 2.2 earthquakes somewhere in the region each day. Therefore, aftershock generation equivalent to an M2.2 in each of the bins in the region is, in some sense, equivalent to a factor of 270 increase in the seismicity, that is, a rate change of 1/0.86 multiplied by the number of bins. Thus dynamic triggering can be a significant feature of the seismic hazard, but it is important to recognize that the power of dynamic triggering comes from its ability to activate a large region and not from a particularly large tendency to triggering in any one spot. Moreover, the dynamic triggering can eventually increase the probability of a large magnitude earthquake on a major fault.
The second observation of a depressed rate of seismicity decay may be the most important result of this study. It is the clearest evidence yet that remote, dynamically triggered events are different from ordinary aftershocks. This measurement also addresses one of the most enduring mysteries of dynamic triggering: Why do the triggered earthquakes persist well after the perturbing seismic waves pass? The observation of low p-values needs to be reconciled with the higher p-values that have been inferred for individual remotely triggered sequences. This reconciliation is helped by realizing that the two different metrics are measuring different quantities. The cumulative number distribution ratios are measuring the seismicity rate directly triggered by the dynamic strain based on the first observed earthquake. The case studies are measuring the decay of earthquake rate following the first local earthquake. The latter likely includes secondary aftershocks of the initially triggered earthquake, which appear to decay quickly. The former is a measure of the triggering process and indicates that the triggered rate decays gradually, which is consistent with previous observations of delayed dynamic triggering that show the multi-day increase in seismicity from background rates (e.g., Parsons et al., 2014). This gradual decay could be indicative of fatigue, local creep, or fluid flow driven by the initial dynamic strain. It cannot be simply a cascade of aftershocks, which would have a more ordinary p-value (Ogata, 1988). Thus, the delayed nature of dynamic triggering includes a secondary process such as creep or fluid flow.
The third observation reinforces the inference that either creep or fluid flow is an important part of dynamic triggering. It seems that dynamic triggering is most prevalent in areas that are prone to creep and/or have known fluid involvement. The tendency of geothermal areas to have triggered behavior is not without exception. Furthermore, the high background seismicity in geothermal areas makes the triggered seismicity more easily observed when triggering occurs as a fixed fraction of the background. As a result, we are more intrigued by the correlation with areas of a creep. The same area showed large triggerability in the previous study of van der Elst and Brodsky (2010) and therefore is one of the most robust dynamic triggering observations. There appears to be a strong relationship between creep and the ability to dynamically trigger.

Conclusions
We successfully reproduced and refined the relationship between the peak dynamic strain and the triggering intensity in Southern California by using new datasets that are independent of previous studies. For the first time, we have established the triggering rate as a function of the directly measured local ground motion. Dynamic strains of 2 × 10 −6 produce triggered events equivalent to the aftershocks of a magnitude 2.2 within each 0.1°×0.1° bin of the study area. Although this is a mild degree of triggering in each local bin, a teleseismic wave can activate the entire region at this level, thus resulting in a significant effect. These empirical relationships can be used to predict seismicity on the basis of ground motions.
The strongest constraint on the physical nature of dynamic triggering comes from the new observation that the triggered seismicity rate decays more slowly than ordinary aftershocks. The slow decay means that the prolonged nature of dynamic triggering is different from ordinary aftershocks. Dynamic triggering is most prevalent in areas of creep and known fluid involvement. Together the observations point toward an important role for creep or fluids in the dynamic triggering process. A cascade of secondary triggering is insufficient to explain the observations. Dynamic triggering includes a distinctive, but still hidden, process.