Rare Occurrences of Non-cascading Foreshock Activity in Southern California.

Earthquakes preceding large events are commonly referred to as foreshocks. They are often considered as precursory phenomena reflecting the nucleation process of the main rupture. Such foreshock sequences may also be explained by cascades of triggered events. Recent advances in earthquake detection motivates a reevaluation of seismicity variations prior to mainshocks. Based on a highly complete earthquake catalog, previous studies suggested that mainshocks in Southern California are often preceded by anomalously elevated seismicity. In this study, we test the same catalog against the Epidemic Type Aftershock Sequence model that accounts for temporal clustering due to earthquake interactions. We find that 10/53 mainshocks are preceded by a significantly elevated seismic activity compared with our model. This shows that anomalous foreshock activity is relatively uncommon when tested against a model of earthquake interactions. Accounting for the recurrence of anomalies over time, only 3/10 mainshocks present a mainshock-specific anomaly with a high predictive power.

Introduction

Large earthquakes are often preceded by an increase in seismic activity, which is then referred to as a foreshock sequence (Bouchon et al., 2013; Dodge et al., 1995, 1996; Jones & Molnar, 1976; Marsan et al., 2014; Reasenberg, 1999). Although these foreshock sequences are often referred to as precursors, a problem is the inherent difficulty to identify earthquakes as foreshocks before the mainshock occurs. In addition, we still do not fully understand the physical mechanisms that generate foreshocks and the reason why they occur. Two competing conceptual models have been proposed (Mignan, 2014). First, a “cascade model” where successive foreshock stress changes contribute to a slow cascade of random failures (possibly mediated by aseismic afterslip) ultimately leading to the mainshock (Ellsworth & Bulut, 2018; Helmstetter & Sornette, 2003; Marzocchi & Zhuang, 2011). Second, a “slow pre‐slip model” where foreshocks are passive tracers of an evolving fault loading process preceding the mainshock rupture (Bouchon et al., 2011; Dodge et al., 1996; Kato et al., 2016). The aseismic versus seismic contributions to the overall moment release during the precursory phase is ultimately what distinguishes these two models. Unfortunately, the aseismic part is generally difficult or merely impossible to estimate from the available observations, and one therefore needs to resort to indirect arguments, often pertaining to the spatial and temporal distribution of the foreshocks. Although recent observations of slow deformation transients lasting days to months before the mainshock favor the triggering of foreshocks by aseismic pre‐slip (Ito et al., 2013; Mavrommatis et al., 2014; Socquet et al., 2017), the aseismic character of such precursory motion is vigorously debated (Bedford et al., 2015; Ruiz et al., 2014). In addition, foreshock sequences are not observed systematically before large earthquakes. However, this lack of systematic precursory observations might partly be due to the incompleteness of current seismicity catalogs (Mignan, 2014; Ross et al., 2019).

The southern California catalog was recently enhanced thanks to the template matching analysis conducted by Ross et al. (2019). The resulting QTM (Quake Template Matching) catalog includes more than 850,000 earthquakes (for the higher choice of threshold, see Section 2.1) in a 10 year‐long period from 2008 to 2017 and is complete down to magnitudes near or below zero for the best resolved regions. Such a high degree of completeness of the QTM catalog motivates the evaluation of the statistical significance of seismic activity preceding large earthquakes in Southern California. By comparing seismic activity before M ≥ 4 earthquakes to a constant and local background rate, Trugman and Ross (2019, T&R from here on) estimated that 72% of mainshocks in the QTM catalog are preceded by a significantly elevated seismic activity. With the same approach using the Southern California Seismic Network (SCSN) catalog, which includes less earthquakes, only 46% of mainshocks were detected with a significantly elevated seismic activity. These results suggest that detailed earthquake detections could bear important information about an impending earthquake. The seismic activity observed in the 20‐day window before M ≥ 4 earthquakes was later re‐evaluated by van den Ende and Ampuero (2020, V&A from here on) to investigate in which cases these increases in seismicity were significant compared to the natural fluctuations of the seismicity rate. In their approach, V&A choose to test seismic activities smoothed at 20 days against a model that accounts for increases in seismicity. In this model, earthquake inter‐event times (IETs) are drawn independently from a gamma distribution. This approach is motivated by the fact that IETs in seismic catalog tends to follow a gamma, rather than an exponential distribution (i.e. T&R's background model) because the gamma distribution is more likely to fit the small IETs observed during clusters of earthquakes. Based on this analysis, V&A estimated that only 33% of mainshocks are preceded in the last 20 days by a significantly elevated seismic activity, coming down to 18% when accounting for temporal fluctuations of such anomalies, i.e., anomalies taking place at random and therefore not specifically related to mainshock occurrences.

For the sake of simplicity, we will now refer to as “foreshock activity” the seismic events observed in the 20 days immediately before M ≥ 4 earthquakes. Although V&A further addressed the significance of elevated foreshock activity in the QTM catalog, we believe that their analysis still underestimates the effect of earthquake clustering. Namely, the random sampling approach of V&A assumes independent IETs, which is an over‐simplification of the actual earthquake clustering observed during individual aftershock sequences. Indeed, during aftershock sequences, IETs are correlated rather than independent. We illustrate this concern in the supporting information (Text S4 and Figure S6) by applying the V&A approach on synthetic ETAS catalogs. In this study, we consider that local earthquake interactions need to be fully accounted for in order to identify foreshock activity that stands out from simple cascades of triggered seismicity.

We extend the studies of T&R and V&A by testing the statistical significance of elevated foreshock seismicity in the QTM catalog, accounting for local earthquake interactions. In this work, we use the temporal Epidemic Type Aftershock Sequences (ETAS) model, in which the seismicity rate at each time is represented by the superposition of a background rate and a rate linked to the aftershock triggering from past events (Ogata, 1988). This model is the simplest that can reproduce both the gamma distribution of IETs (Saichev & Sornette, 2007) and their correlation during aftershock sequences. After selecting mainshocks using criteria similar to T&R and V&A, we extract ETAS parameters from the QTM catalog in the vicinity of each mainshock. We then compare the foreshock activity with ETAS predictions accounting for past seismicity. We find that the number of instances of anomalously elevated foreshock seismicity is significantly reduced when accounting for earthquake interactions (about 19% compared to 33% and 72%, respectively, in V&A and T&R). Moreover, out of these 10 cases, only 3 appear to be exclusively related to the subsequent occurrence of the mainshock.

Data and Methods

Mainshock Selection

We noticed that the full QTM catalog used by T&R and V&A suffers from episodic bursts of false detections, that occur due to too low a detection threshold (threshold fixed at 9.5 times the median absolute deviation (MAD) of the stacked correlation function). These bursts are easy to identify as they start or end at midnight, which is due to the MAD computation being performed over 24 h long period starting at 00h00 UTC. To avoid any contamination of our analysis by such artifacts, we instead use the higher quality QTM catalog with a detection threshold at 12 times the MAD, for which these transients vanish or are strongly attenuated. In order to provide a fair comparison with previous results, we also present our analysis performed on the full catalog in the supporting information (Text S5 and Figures S7 and S8).

Using the higher quality QTM catalog, we then extract our own set of mainshocks with selection criteria similar to those used in T&R: A mainshock must have magnitude M ≥ 4, and must occur from 2009/01/01 to 2016/12/31 within the geographic coordinates ranges [32.68°N, 36.2°N] and [118.80°W, 115.4°W]. To be selected, a mainshock must be preceded by at least 10 earthquakes with no larger magnitude event in the year before and within a 20 × 20 km2 horizontal box around its epicenter. 53 earthquakes were selected as mainshock according to these criteria. For each selected mainshock, we extract a 10‐year long local catalog that includes all the seismicity observed within the 20 × 20 km2 box with no depth cutoff.

We evaluate for each local catalog the local magnitude of completeness M and remove all events with a magnitude M < M . We must acknowledge that removing all earthquakes of the QTM catalog below M may remove potentially interesting features, but we consider that such features cannot be properly interpreted because they might reflect variation of the detection capability of the network and not real fluctuations of the seismicity rate. Therefore, to achieve a trade‐off between completeness and retaining as many earthquakes as possible, we estimated manually the local M as either the maximum of the local Gutenberg‐Richter (G‐R) frequency‐magnitude distribution if this distribution decays smoothly for larger magnitudes, or the magnitude at which a notable break in slope is observed. Figure S1 shows the 53 local Gutenberg‐Richter frequency‐magnitude distributions and the corresponding estimated M values.

Inversion of ETAS Parameters

The ETAS model has two main ingredients: First, a background term which is time‐independent and follows a Poisson process; second, a triggered term that depends on the past earthquake activity. The conditional intensity of the ETAS model (Ogata, 1988; Zhuang et al., 2012) is:

where μ is the time‐independent background seismicity rate. The sum in the right hand side of Equation  1 describes the expected aftershock seismicity rate at time t triggered by all previous events. A and α are constant parameters describing respectively the global aftershock productivity of the region and the magnitude dependence in the number of triggered events. M is the magnitude of completeness whereas c and p are the parameters of the Omori‐Utsu law describing the time‐decay in the aftershock seismicity rate. Therefore, in ETAS‐like catalogs, temporally clustered seismicity only emerges from cascades of aftershocks.

For local catalogs associated with each mainshock, we fit the temporal ETAS model by maximizing a likelihood function with an Expectation ‐ Maximization (EM) algorithm (Veen & Schoenberg, 2008). We estimate parameters A, c, p, α, and μ in equation (1) (all parameter values can be found in the supporting information). We run a first inversion where the ETAS parameters are constrained to be positive. We note that most α values are close to one. Larger α values are actually expected according to window‐based methods (Felzer et al., 2004; Helmstetter, 2005), as well as following the argument that Bath's law, i.e., the fact that the difference in magnitude between the mainshock and its largest aftershock is independent of the mainshock's magnitude, requires that α = β = b ln 10 (Davidsen & Baiesi, 2016 and references therein). Moreover, it has been shown that α estimates are particularly prone to model errors (e.g., Hainzl et al., 2008, 2013) and censoring effects (Seif et al., 2017; Sornette & Werner, 2005). Nandan et al. (2017) found that the α value is expected to vary between 1.7 and 2.2 when considering a larger portion of California and a longer period than the QTM catalog. A α value close to 2 may thus represent a more realistic value of the aftershock productivity for Californian earthquakes. Therefore, we perform a second inversion where we impose that α = 2. We thus obtain two sets of ETAS parameters (referred to as “α free” and “α = 2” sets) to model the seismicity of local catalogs around each mainshocks. We also evaluate in the supporting information the sensitivity of our results to the uncertainty in ETAS estimates for both sets of parameters (cf., Text S3 and Figures S4–S5).

Detection of Seismicity Anomalies Based on the ETAS Model

We test the null hypothesis H 0 that the number of events observed in 20 days is smaller than or equal to the number of events predicted by the ETAS model for both sets of parameter estimates. If H 0 is rejected for both estimates, we assume that an anomalously high seismicity is detected in the window, suggesting that a mechanism other than simple ETAS cascading is required to explain the 20 day earthquake activity. The conditional intensity function in equation (1) allows to directly compute an expected seismicity rate at any time t from the set of ETAS parameters (A, c, p, α, and μ) and the knowledge of past seismicity (t  < t, M ). By integrating this modeled seismicity rate, we can compute the expected number of earthquakes in a time interval T:

Here, we set T = 20 days similar to T&R, which choice was also adopted by V&A. We compute over 20 day sliding windows, with a 1 day shift between two consecutive windows, and covering the full time range of the QTM catalog (i.e., 10 years). For all local catalogs around each mainshock, we then obtain two time‐series of generated using the two sets of inverted ETAS parameters (α free and α = 2). Knowing , the probability of actually observing N obs earthquakes in a given 20 day time‐interval is given by the Poisson distribution with mean :

We then define the probability of observing at least N obs events over 20 days for the null hypothesis as:

Following T&R and V&A, we use the probability threshold p < 0.01 to reject the hypothesis H 0 that N obs is in agreement with the expected number of events . A small p‐value would therefore correspond to anomalously elevated seismicity rate compared with ETAS predictions.

Results

The detection of seismicity rate anomalies in a 20 days sliding window is illustrated in Figure 1 for the seismicity located in the vicinity of four mainshocks. For each mainshock, the top subplot shows the time‐evolution of p‐values measured for the two sets of ETAS parameters (α free and α = 2) while the bottom subplot shows the observed seismicity (i.e., magnitude vs. time). For the two examples on top (Mainshock IDs 10832573 and 37301704), we notice that the 20 days foreshock activity is consistent with ETAS predictions with a p‐value above 0.01 in the last 20 days window prior to the mainshock. In these cases, our null hypothesis H 0 cannot be rejected with a confidence of 99%. The two examples on the bottom (Mainshock IDs 14898996 and 37299263) show p‐values that are below 0.01 before the mainshock for both ETAS estimates. In these cases, the observed foreshock seismicity is higher than the expected ETAS cascading seismicity with a confidence level of at least 99%.

Figure 1

The 20‐day sliding window analysis for four examples of mainshocks (black star at t = 0) and their local catalog. Mainshocks IDs are (a) 10832573, (b) 37301704, (c) 14898996 and (d) 37299263 (Top graphs) probability p that ETAS explains the observed seismicity, computed for the two sets of ETAS parameters α = 2 and α free. The p‐value for the last 20‐day window prior to the mainshock is shown with a thick square. The significance threshold of p = 0.01 is shown with the horizontal dotted line (Bottom graphs) magnitude versus time for the local catalogs in the 20 × 20 km2 box around each mainshock. The right inset is a zoom around the foreshock window. ETAS, Epidemic Type Aftershock Sequences.

In total, we find that 10 out of 53 mainshocks are preceded by an anomalously high 20 days activity with respect to ETAS predictions. Therefore, these mainshocks are likely preceded by complementary aseismic processes other than cascades of aftershocks. However, this result must be taken in perspective with the overall ability of the ETAS models to explain fluctuations in seismicity rates over the entire catalog. As pointed out by V&A, the predictive power of an anomalously high foreshock activity is reduced if seismicity anomalies are frequently detected without being followed by a large event. The significance of an anomalously high foreshock activity being predictive of future large events should therefore be assessed given the overall ability of ETAS predictions to explain the seismicity in the vicinity of the mainshock. For example, in the case of mainshock ID 14898996 in Figure 1c, ETAS predictions are unable to explain the observed seismicity at several occasions during the course of the catalog. Our null hypothesis H 0 is thus rejected for numerous 20 days windows with p‐values smaller than the p‐value of the foreshock window. On the other hand, Figure 1d shows that mainshock 37299263 presents an anomalously high seismicity rate almost exclusively in the 20 days preceding the mainshock. Such an elevated seismicity rate is thus highly correlated with the mainshock occurrence. We believe that the uniqueness of the anomaly observed before mainshock ID 37299263 is more likely to evidence predictive non‐cascading mechanisms than mainshock ID 14898996.

Therefore, to quantify the significance of detected foreshock anomalies, we compare p‐values in the foreshock window with the distribution of p‐values over the entire 10‐year catalog. For each mainshock, an anomalous foreshock activity is considered mainshock‐specific if , the proportion of 10‐year p‐values lower or equal than the foreshock p‐value, is less than 1%. This arbitrary threshold of 1% allows to discriminate between catalogs with frequent anomalous activities and those with foreshock activities that correspond to the strongest anomalies of their region. This is summarized in Figure 2b. Using such temporal specificity criterion, we identify that 7 out of the 10 anomalous foreshock activity already mentioned occur in regions with recurrent seismicity anomalies stronger than the foreshock one. Therefore, we argue that only 3 out of 53 mainshocks present a clear mainshock‐specific anomalous activity. We note that this final selection is highly dependent on the choice of the threshold. Figure 2b shows that all 10 selected sequences present less than 10% of 20‐day windows over 10 years below the foreshock window p‐value. The final selection of 3 out of 53 mainshock is therefore more like a refined selection of mainshocks with a local seismicity that best fit ETAS with a notable exception during foreshock time ranges.

Figure 2

(a) Receiver Operating Characteristic (ROC) curves for our detection of anomalous foreshock windows. For a varying threshold p‐value p thresh, the curves show the proportion of foreshock windows below p thresh against the proportion of non‐foreshock windows below p thresh. ROC curves are drawn for both the full set of 53 local catalogs and the set of 43 catalogs left after removing the 10 anomalous sequences of section 3 (with p < 0.01). We also include the ROC curve corresponding to the average of 53 sets of 1000 ETAS simulations computed using the α free ETAS parameters obtained in Section 2.2. Note that ETAS simulations display a curved ROC, the departure from the “no‐gain” line being particularly clear when considering large p thresh values. This departure is weak for p thresh ≤ 0.01, with a gain of about 2 at maximum (p thresh = 0.01). (b) Proportion of windows with a p‐value lower or equal to the 20‐day foreshock window p‐value, among all 20‐day windows over 10 years. The proportion is shown here for the 10 anomalously high foreshock activity and for the two ETAS estimates. We consider an anomalously high foreshock activity to be specifically related to a mainshock if is below 0.01 for both ETAS estimates. Here, we identify three foreshock anomalies that are specific to subsequent mainshocks for both sets of ETAS parameters. Note that is significantly sensitive to the value of α. Labels preceded by a star are mainshock IDs of the two anomalously high foreshock activity detected with the declustering approach. ETAS, Epidemic Type Aftershock Sequences.

We complement this analysis by declustering the local catalogs. The probability ω that earthquake i is a background earthquake is defined as , and can be calculated once the ETAS parameters are estimated. We then simply count the numbers of background earthquakes as the sums of ω in 20 day long windows. We denote N 0 this count for the last 20 days prior to the mainshock, and by N all the counts for all the time windows before the mainshock (not just the last one). Following the same rationale that stimulated our previous analysis, we first compare N 0 to the Poisson distribution with a mean equal to the mean of N, select the mainshocks for which for the two sets of ETAS parameters (first test), and finally check whether these selected sequences display other anomalously strong bursts of background earthquakes by computing the probability that N can be greater than N 0 (second test). We finally select those short‐listed mainshocks for which the latter probability is less than 0.01 (again, for the two sets of ETAS parameters). Figure 3 shows the results of this declustering approach. Only mainshocks 14598228 and 14600292 are preceded by an anomalously high foreshock activity (first test) according to this declustering approach. According to our second test, these two anomalies are also specific to the subsequent mainshock occurrences (i.e., p‐value ≤ 0.01). These two foreshock sequences were also identified in our previous approach based on the predicted number of events according to the ETAS model. The difference in results between the declustering approach and the former method is due to the fact that declustering only leaves a small number of background earthquakes, and therefore has a strong tendency to significantly lower the p‐values.

Figure 3

(a,b,c) Number of earthquakes in 20 day long windows counting (top) all earthquakes and (bottom) background earthquakes only, for three selected mainshocks. The number for the last window prior to the mainshock is shown with a thick square. The dashed lines show, for the two sets of ETAS parameters (free α in red, α = 2 in blue) the limit over which the Poisson probability becomes less than 0.01. (d) Probability P(N > N 0) that the last 20 days are anomalously active compared to the past, for the two sets of ETAS parameters; the sequence is selected as a mainshock‐specific anomalous activity after declustering if this probability is less than 0.01 (second test) and if N 0 is above the dashed line (first test). Mainshocks 14598228 and 14600292 correspond to indices 0 and 1 on this graph, and are the only mainshocks with both probabilities less than 0.01. All indices can be linked with their mainshock ID thanks to Table S2.

Discussion

We use the highly complete QTM catalog of Ross et al. (2019) for Southern California to further investigate the significance of anomalous high foreshock activity previously reported by T&R and V&A. As mentioned before, those studies did not fully address whether the temporal clustering of earthquakes observed during aftershock sequences is a possible explanation for the observed elevated foreshock activities. This clustering is considered as one of the possible origins of the high seismic activity observed before large earthquakes (Ellsworth & Bulut, 2018; Helmstetter & Sornette, 2003; Marzocchi & Zhuang, 2011). In practice, small M < 4 earthquakes trigger small aftershock sequences during which a larger M > 4 event is more likely to occur than at more quiet times. In this regard, high activity preceding a mainshock can naturally stem from such earthquake interactions and cascading without necessarily requiring an external pre‐slip phenomenon. To address this concern, we use the ETAS model to discriminate which instances of QTM foreshock activities exhibit higher seismicity rates than expected from earthquake interactions.

We first assess the probability p that a given 20 day foreshock sequence can be explained by ETAS earthquake clustering. Using p < 0.01 as a threshold, our results indicate that ∼19% (10 out of 53) of mainshocks are preceded by increases in seismicity higher than 99% of the earthquake rates predicted by ETAS. The 20‐day temporal evolution of these 10 anomalous foreshock sequences is detailed in Text S2 and Figure S9. In a second step, we further distinguish 3 out these 10 cases as being specific to the subsequent mainshock, i.e., the chance to see such a significant increase of activity occurring at random is less than 1%. The anomalously high seismicity of these 3 foreshock sequences is thus highly correlated with the M ≥ 4 mainshock occurrences and likely to be controlled by aseismic nucleation processes. We notice that this number (3 out of 10) would raise to 5 if accepting a threshold at 1.5% rather than 1%, cf. Figure 2b. The complementary declustering approach restricts the anomalously high foreshock activity to only two mainshock‐specific sequences. A possible over‐estimation of the background rate can be a cause for this more conservative selection. Even if the definitions of an anomalously elevated seismicity differ, Mainshock IDs related to the anomalously high foreshock activities detected in T&R, V&A and this study can be found in Table S1. The Southern Californian location of these sequences are also compared in Figure S10.

We must emphasize that these results, along with those of T&R and V&A, likely depend on the initial choice of focusing on foreshocks in a 20 day period prior to each mainshock. Using a longer or shorter time‐window may therefore provide different results. Moreover, the fixed 20 × 20 km2 horizontal spatial window used in this study implies that all events in this box are evaluated with the same weight. This can artificially enhance the triggering role of foreshocks that are relatively far from the mainshock. The ETAS model used here would need to be extended to a space‐time model in order to exploit the distance between earthquakes and to help to discriminate such cases (Zhuang et al., 2011, for a review). While this development does not appear over complicated, and was already investigated in Seif et al. (2019), the addition of several model parameters and the use of an isotropic spatial kernel for which no clear consensus exists (Moradpour et al., 2014) is likely to undermine the robustness and significance of the resuts.

The exact number of detected foreshock anomalies obviously depends on the significance threshold that we have fixed to p < 0.01 following T&R and V&A. To assess the impact of this arbitrary choice, we evaluate how the proportion of detected anomalous high foreshock activity changes as a function of the p‐value threshold p thresh. This result is compared with the proportion of windows that have p < p thresh without being followed by a mainshock (i.e., false positives). We thus compute the Receiver Operating Characteristic (ROC) curve as shown in Figure 2a. If the occurrence of anomalously elevated activity was not a sign of an incoming mainshock, then the ROC curve would follow a 1 to 1 straight line (hereafter referred to as the no‐gain line). We find that there is positive correlation between preceding high activity and mainshock occurrence: the information gain is measured by the ratio of true positives over false positives, which is practically constant and close to 6 for p thresh ≤ 0.05. We however notice that significant departure from this no‐gain line also exists in ETAS simulations computed with the same 53 sets of parameters as obtained for the local catalogs. Figure 2 shows that a large p thresh (i.e., p thresh > 0.01) allows to detect anomalous foreshock activities (i.e., a positive gain) in ETAS simulations, even though there is by definition no pre‐slip in this model. This is caused by the clustering properties of the model: In the rare occasions where the observed number of earthquakes N obs in a window largely exceeds the expected number , then the occurrence of earthquakes immediately after this window is more likely, including the occurrence of a mainshock. As an effect, the ROC curve departs from the no‐gain line. We however notice that there is no information gain on the magnitude of the forthcoming earthquakes, as expected. We conclude that choosing too large a value of p thresh may lead to the detection of “foreshock cascades” prior to mainshocks, which are not related to aseismic processes (e.g., pre‐slip). According to our simulations, p thresh = 0.01 appears as an acceptable threshold to discriminate a cascading‐like seismicity from other processes that would also enhance the seismic activity: at p thresh = 0.01, the information gain for ETAS is about 2, compared to about 6 for the observed seismicity (cf., p thresh = 0.01 in Figure 2). This additional gain is mostly controlled by the 10 sequences we found to be anomalous: Quite obviously, removing them from the calculations implies that the ROC curve is equal to zero at p thresh = 0.01. Therefore, these 10 anomalous foreshock sequences suggest the existence of a precursory pattern before some M ≥ 4 earthquakes stronger than expected from ETAS simulations.

Our results strengthen previous reports that earthquake activity precursory to mainshocks can sometimes deviate from simple clustering properties (as modeled by ETAS; Lippiello et al., 2019; Seif et al., 2019). Our approach is however different. For example, compared to Seif et al. (2019), we seek to explain the last 20 days prior to mainshocks knowing all past seismicity (including activity in the last 20 days), by comparing what number of earthquakes would be “normally” expected (in the sense of ETAS) to the observed number. In contrast, Seif et al. (2019) compared observations to the number of foreshocks predicted by ETAS simulations not constrained by past seismicity. Our method is indeed close to the residual analyses of Ogata (1988), (1989), (1992), and (2003), which is here performed individually on a set of 53 mainshocks thanks to the improved completeness of the QTM dataset.

Conclusions

According to our analyses, the low magnitude of completeness of the QTM catalog does not warrant the detection of aseismically driven foreshock sequences in the 20 days window preceding isolated mainshocks. More than 80% of mainshocks are preceded in the last 20 days by activity exhibiting seismicity rates that are consistent with ETAS predicted rates, even when the magnitude of completeness is as low as M  = 0. For these cases, earthquake interactions and local stress changes are a good candidate to explain all observed increases in seismicity rates prior to the mainshock. We find 10 mainshocks that are preceded in the last 20 days by a significantly high seismic activity. These cases show seismic activity that significantly differ from ETAS cascades, and are thus likely controlled by aseismic processes. Among those 10 cases, we distinguish three cases that exhibit non‐ETAS like seismicity that is very likely specifically related to the mainshock; these three cases are the best evidences of a possible nucleation phase.

High quality earthquake datasets complete to low magnitudes are in any case required to pursue and develop efforts for understanding when and where aseismic pre‐slip can lead to a large shock. Foreshocks remain the best observable to study preparatory processes, if they exist (Nakatani, 2020). First, increasing the location accuracy and the number of small earthquakes substantially improves the statistical significance of any test conducted to assess the reality of pre‐slip processes, when comparing to the cascade (null) hypothesis. Second, the availability of large dataset. allows to increase the number of potential mainshocks to be analyzed, hence offering more robust conclusions. Finally, we suggest that pre‐slip seismicity analysis should be evaluated along other near‐fault observables (such as GPS data (Socquet et al., 2017), strainmeter data (Roeloffs, 2006), variations in groundwater level or flow rate (Roeloffs, 1988), radon emission rate (Ghosh et al., 2009), changes in seismic velocities as imaged by pairwise seismic station cross‐correlation functions (von Seggern & Anderson, 2017)) whenever available, to independently assess any possible aseismic mechanisms at work during the preparation of large earthquakes.

Supporting Information S1

Click here for additional data file.

Table S1

Click here for additional data file.

Table S2

Click here for additional data file.