1 Introduction

The existence of precursory signals before large earthquakes remains a central question in seismology (Cicerone et al., 2009; Mignan, 2014; Peng & Lei, 2025; Pritchard et al., 2020; Scholz et al., 1973; Wyss, 1991). Foreshocks, in particular, have been widely investigated as potential precursors, yet their prognostic value is still controversial (Felzer et al., 2004; Petrillo & Lippiello, 2020, 2023; Seif et al., 2019; Zaccagnino et al., 2024). While laboratory experiments and numerical models often reveal accelerating slip prior to failure (Dieterich & Kilgore, 1996; Kaneko et al., 2016; Rubin & Ampuero, 2005), translating these findings to natural faults remains challenging (Hulbert et al., 2019; Passelègue et al., 2017). One key limitation is that laboratory studies typically assume homogeneous fault conditions, whereas natural faults exhibit strong structural and frictional heterogeneity—including localized seismic and aseismic slip, variations in frictional properties, and cascading failure processes—all of which can critically influence rupture nucleation and precursor expression.

Recent advances in seismic and geodetic monitoring have enabled the identification of transient fault activity before some large earthquakes, including slow slip events and foreshock sequences (Beaucé et al., 2023; E. E. Brodsky & Lay, 2014a; A. Kato & Ben-Zion, 2021; Martínez-Garzón & Poli, 2024). Notably, a recent study by Bletery and Nocquet (2023, 2025) suggested the existence of systematic precursory slip in geodetic data, sparking both interest and debate in the scientific community (Bradley & Hubbard, 2023, 2024). However, such transient signals are not systematically observed, and similar activity has often been detected without an ensuing mainshock (Behr & Bürgmann, 2021; Bletery & Nocquet, 2020; Gomberg et al., 2010; Obara & Kato, 2016). This lack of consistency makes it difficult to determine whether seismic sequences exhibit diagnostic signals—observable in seismic or geodetic data—that indicate imminent fault failure.

The complexity of fault behavior is often associated with spatial variations in frictional properties, with faults containing regions that either strengthen or weaken as slip rates change (Avouac, 2015). During interseismic periods, rate-strengthening fault patches creep steadily under tectonic loading, while rate-weakening patches remain locked and accumulate stress (N. Kato, 2020; Kazemian et al., 2015; Pelletier, 2000; Lippiello, Petrillo, Landes, & Rosso, 2019; Lippiello et al., 2021; Petrillo et al., 2020). When these locked patches experience stress levels that exceed their local strength, they may undergo seismic slip, perturbing adjacent regions and potentially triggering cascading failure along the fault (McLaskey, 2019; Schwartz & Rokosky, 2007). Large ruptures are generally thought to nucleate within critically stressed fault patches; their growth into major earthquakes depends on the rupture's ability to propagate through heterogeneous fault regions and activate adjacent areas, ultimately leading to collective failure. While rate-and-state friction laws provide a widely used framework for modeling fault slip, they offer only a simplified description of natural rupture processes. Near failure, rupture dynamics are increasingly governed by the structural complexity and topological organization of a fault, which influence how stress is transferred and redistributed across the fault. In this regime, dynamic fracture and long-range interactions can dominate over local frictional properties, suggesting that rupture is better understood as an emergent collective phenomenon rather than the sum of isolated patch instabilities (Fineberg & Marder, 1999; Lee et al., 2024; Reches & Fineberg, 2023). Since direct stress measurements on faults are not feasible, indirect methods are necessary to assess how close a fault is to failure. Geodetic observations provide constraints on slip deficit, but their spatial resolution is often insufficient to resolve fine-scale variations in fault stress (Cohen, 1999; Okada, 1985). Consequently, seismic data remain a critical tool for investigating earthquake nucleation and rupture dynamics.

Seismic waveforms contain rich information about asperities—localized regions where stress accumulates and is released (Lay & Kanamori, 1981; Lei, 2003). Previous studies have shown that the frequency content of seismic waves provides insight into asperity size: smaller asperities radiate higher-frequency waves due to their rapid energy release, whereas larger asperities release energy more slowly, producing lower-frequency signals (Aki, 1967; Kumagai et al., 2012; Sato & Hirasawa, 1973). However, the potential of waveform analysis for identifying evolving fault conditions remains underexplored. In this study, we focus on foreshock waveforms, analyzing 68 M6+ earthquakes in selected regions since 2011 that occurred near and shortly before a larger event (Figure 1). Our results show that foreshocks preceding large earthquakes often exhibit characteristic waveform patterns—marked by rapid fluctuations and elevated energy release—that differ from isolated events of similar magnitude. These observations suggest that waveform analysis, particularly through the lens of early envelope behavior, may offer indirect constraints on the state of stress of faults and their potential to host a larger rupture. This has direct implications for the ongoing debate on whether foreshocks reflect preparatory processes or stochastic clustering.

Although foreshocks are frequently observed before large earthquakes, they also occur prior to smaller events, making their predictive significance uncertain (E. E. Brodsky, 2011; de Arcangelis et al., 2016; Lippiello et al., 2012, 2017; Moutote et al., 2021; Ogata & Katsura, 2014; Petrillo & Lippiello, 2020, 2023; Seif et al., 2018; Shearer, 2012; Trugman & Ross, 2019; van den Ende & Ampuero, 2020). Two main models have been proposed. The cascade model (Fukao & Furumoto, 1985; Helmstetter & Sornette, 2003; Marsan & Enescu, 2012) views foreshocks as typical earthquakes that can occasionally trigger a larger rupture through stress transfer, without reflecting any special physical state. In contrast, the pre-slip model (Dodge et al., 1996; Ellsworth & Beroza, 1995) suggests that foreshocks mark the onset of a preparatory process, where localized aseismic slip gradually loads adjacent fault patches until dynamic rupture occurs. Support for this model includes geodetic observations of slow slip between foreshock and mainshock (Ito et al., 2013; A. Kato et al., 2012; E. E. Brodsky & Lay, 2014b), and seismicity changes such as a decreasing Gutenberg-Richter b-value prior to large events—opposite to the increase typically seen after mainshocks (Gulia & Wiemer, 2019; Gulia et al., 2018). Although the interpretation of b-value trends remains debated (DeSalvio & Rudolph, 2021), this pattern has been documented in multiple foreshock–mainshock sequences (Chan et al., 2024). For further context, see recent reviews (Martínez-Garzón & Poli, 2024; Peng & Lei, 2025).

In this study, we introduce a novel approach to foreshock analysis by identifying anomalous patterns in waveform envelopes. We show that these anomalies, quantified by our newly developed Q index, systematically appear in foreshocks preceding large earthquakes. Our results provide insight into the physical processes governing foreshock sequences and contribute to the broader discussion on the seismic signals that may precede major ruptures.

2 Methods and Results

The central quantity in our analysis is the envelope function μ(t), defined as the logarithm of the envelope of the ground velocity (Kanasewich, 1981) (see Supporting Information S1) (Figure 2). Prior to an earthquake, μ(t) fluctuates around a background level μ_B, before rising sharply and reaching a maximum value μ_M at time t_M, corresponding to the perceived magnitude at the recording station (Lippiello et al., 2016).

The post-mainshock evolution of μ(t) is strongly influenced by aftershock activity (Lippiello et al., 2016; Lippiello, Cirillo, et al., 2019; Lippiello, Petrillo, Godano, et al., 2019; Peng et al., 2007; Sawazaki & Enescu, 2014). Figures 2a and 2b illustrate two contrasting cases: the 27 December 2019 M_w6.1 Kissamos earthquake, which was followed by no significant aftershocks (μ ≤ μ_B), and the 20 July 2017 M_w6.6 Kos earthquake, which exhibited multiple large-amplitude secondary peaks consistent with elevated aftershock productivity in the hours following the mainshock. Following its peak at t = t_M, the envelope function μ(t) of the M_w6.1 Kissamos earthquake decays rapidly and stabilizes at the background level, μ(t) ≃ μ_B, consistent with an absence of aftershock activity. In contrast, for the M_w6.6 Kos earthquake, μ(t) does not decay directly to μ_B, but instead remains elevated, fluctuating around a plateau level μ_L ≫ μ_B. The background level is only recovered after several weeks (Lippiello et al., 2016; Lippiello, Petrillo, Godano, et al., 2019; Shearer & Senobari, 2025), well beyond the time window shown in Figure 2b. This plateau is a hallmark of intense aftershock activity (Lippiello et al., 2016): when an aftershock occurs at time t₁ > t_M with perceived magnitude μ₁, the envelope remains elevated at μ(t₁) = μ₁, preventing an immediate return to μ_B.

If a second aftershock with perceived magnitude μ₂ occurs at t₂ > t₁, the envelope rises again, reaching μ(t₂) = μ₂ > μ_B. When aftershock productivity is high, the interval (t₂ − t₁) between successive aftershocks is short, preventing μ(t) from dropping below the plateau level μ_L ≃ min(μ₁, μ₂) > μ_B. Numerical simulations (Lippiello et al., 2016) show that the number of aftershocks increases exponentially with μ_L − μ_B. More specifically, the expected number of aftershocks (μ_i > μ_M − 2) following a mainshock of perceived magnitude μ_M scales as n_aft ∝ 10^−δμ, where δμ = (μ_M − μ_B) − 2(μ_L − μ_B). A detailed estimate of aftershock probability from μ(t) in the first minutes after t_M is provided in Lippiello, Petrillo, Godano, et al. (2019).

The envelope function after the M_w9.1 Tohoku earthquake of 11 March 2011 (Figures 2c and 2d) exhibits the expected high plateau level μ_L, reflecting intense aftershock activity. However, the behavior observed after the M_w7.3 Tohoku foreshock (Figures 2c and 2d) is strikingly different. In the first 8 minutes, μ(t) returns to values close to μ_B, resembling the no-aftershock pattern seen in the M_w6.1 Kissamos earthquake (Figure 2b). Unlike Kissamos, however, the Tohoku foreshock envelope undergoes abrupt rises and falls, forming an anomalous sawtooth pattern. Large peaks in μ(t) correspond to significant foreshocks (μ_i ≃ μ_B + 4), whereas valleys with μ(t) ≃ μ_B indicate periods with very few μ_i > μ_B events. This contradicts the Gutenberg-Richter law, which predicts thousands of events with μ_i ≃ μ_B + 1 for every event with μ_i ≃ μ_B + 4. It also challenges the cascade model, which interprets fore-mainshock sequences as typical triggering processes. In this model, a mainshock is simply an earthquake that happens to be larger than its triggering event, that is, the foreshock. The occurrence of an aftershock larger than the triggering event is expected in periods of intense aftershock activity (Fukao & Furumoto, 1985; Helmstetter & Sornette, 2003; Spassiani et al., 2024). Consequently, the cascade model would predict an envelope function μ(t) with a high plateau level μ_L, reflecting a large aftershock count n_aft. However, the observed envelope profile (black curve in Figure 2d) instead reveals the opposite trend: the Tohoku foreshock exhibited a limited ability to trigger small aftershocks.

This contrast between foreshock and mainshock envelope behavior is consistently observed across different seismic stations located at varying azimuths and distances from the rupture (see Text S6 in Supporting Information S1). Examining the evolution of μ(t) in other foreshock-mainshock sequences reveals a similar pattern. Figures 2f and 2h show the envelope function μ(t) for the M_w8.1 1 April 2014 Iquique earthquake and the M_w7.1 6 July 2019 Ridgecrest earthquake, respectively. In both cases, the mainshock envelope exhibits a high plateau level, μ ≫ μ_B, indicating intense aftershock activity. In contrast, the envelopes of their large foreshocks (Figures 2f and 2h) display the sawtooth pattern similar to that observed in the Tohoku foreshock. This recurring behavior across multiple earthquake sequences suggests that waveform envelopes may provide valuable insights into foreshock-mainshock dynamics.

2.1 Quantifying Foreshock Waveform Anomalies

To systematically evaluate the observed differences in the envelope function (Figure 2), we introduce a quantitative measure to assess the degree of anomaly in a given envelope. This measure enables us to distinguish between normal envelopes and those exhibiting a sawtooth-like pattern.

Under normal conditions, the peak distribution P(μ) follows the Ishimoto-Iida law (Ishimoto & Iida, 1939), expressed as $$P\left({\mu }_{i}\right)\sim 1{0}^{-\beta {\mu }_{i}}$$, where β is approximately equal to the b-value of the Gutenberg-Richter law. Consequently, the probability of an earthquake with perceived magnitude μ_M being followed by a larger event is proportional to 10^−βμ. However, in foreshock sequences, the number of observed earthquakes n_obs with high peak values (μ_i > μ_M − 2) may significantly exceed the number of expected aftershocks estimated assuming that no further fault segments or systems become unstable as a consequence of the stress redistribution promoted by previous events n_aft ∝ 10^−δμ, which is estimated from the plateau level μ_L. The ratio n_obs/n_aft should therefore exhibit anomalously large values during foreshock sequences. To quantify this, we define the anomaly index $$Q=\left({n}_{\mathit{obs}}/{n}_{\mathit{aft}}\right)1{0}^{-\beta {\mu }_{M}}$$ (see Text S2 in Supporting Information S1). The number of observed peaks, n_obs, is determined following the method by Peng et al. (2007) for μ(t) in the first 45 min after the mainshock. Given that the number of significant peaks is typically small and well separated, a manual count ensures accurate identification of local maxima μ_i exceeding the threshold μ_i > μ_M − 2.

To test this approach systematically, we analyze all M_w6+ earthquakes worldwide since 2011, including additional fore-mainshock pairs from Western South America (Figure 1a), and all M_w6+ events since 2015 in the three primary study regions shown in Figures 1b–1d. These are defined as an M_w6+ event followed within 10 days and 50 km by an event of equal or greater magnitude, with additional quality criteria applied (Text S3 in Supporting Information S1). The complete list is provided in Text S5 in Supporting Information S1. Most events are located in regions with dense seismic networks, such as inland Japan, Central Europe, and Western North America, with additional cases from Central and South America. The data set includes 68 M6+ earthquakes. Of these, 11 were followed within 10 days by a larger earthquake at an epicentral distance of less than 50 km, classifying them as foreshocks. The full selection and classification procedure is described in detail in the Text S5 in Supporting Information S1 together with the seismic network and station chosen for the download. We initially assume β = 1 in the definition of Q and classify earthquakes as anomalous if Q ≥ Q_th and regular if Q < Q_th. The value of Q for all 68 events is provided in Text S5 in Supporting Information S1 and plotted in Figure 3a. Setting Q_th = 0.18, we find that 10 of the 11 foreshocks in the data set are correctly classified as anomalous (Q ≥ Q_th), while 53 of the 57 non-foreshock events are correctly classified as normal (Q < Q_th). Thus, only one foreshock is a missed alarm, and seven of 64 normal events are false alarms. Figure 3 demonstrates that anomalous earthquakes are strongly associated with foreshocks.

To further validate this classification, we use a receiver operating characteristic (ROC) diagram. This approach allows for a variable threshold Q_th instead of a fixed one, enabling us to assess the accuracy of our classification. Specifically, we apply a binary decision rule, triggering an alarm when Q ≥ Q′, and compute the true positive rate (TPR), defined as the fraction of detected foreshocks among all foreshocks in the data set. We also evaluate the false alarm rate (FAR), defined as the fraction of non-foreshock events incorrectly classified as foreshocks. By varying Q′, we obtain the ROC curve (Figure 3b), which shows points near the ideal prediction scenario (TPR = 1, FAR = 0) and significantly above the random prediction line (TPR = FAR). The null hypothesis of random classification is rejected with confidence level cl > 99%, supporting the conclusion that anomalous earthquakes are statistically associated with foreshocks. Further differences between foreshocks and mainshocks emerge in the peak distribution P(μ) (Text S4 in Supporting Information S1). For instance, in the Tohoku sequence, the β-value after the M_w7.3 foreshock (β = 0.4 ± 0.1) is significantly lower than after the M_w9.1 mainshock (β = 1.0 ± 0.1). This change in β is consistent with the b-value decrease observed by Gulia and Wiemer (2019) and is observed in other foreshock-mainshock sequences (Figure 2, upper panels).

Incorporating β from P(μ) into the Q index, rather than assuming β = 1 (as in Figure 3), could improve foreshock identification. However, unlike parameters such as μ_L, μ_B, and μ_M, which can be estimated within the first hour of the envelope function, determining β requires a longer temporal window and a sufficient number of peaks for a reliable fit using the Ishimoto-Iida law. When β is available, the Q index improves foreshock discrimination, correctly identifying all foreshocks as anomalous with only three false alarms (Text S8 in Supporting Information S1). The β values used in this study were estimated using the maximum likelihood method introduced by Aki (1965). Although more advanced estimators, such as the b-more-positive method (Lippiello & Petrillo, 2024), can yield more accurate values under certain conditions the aim of our study is not to optimize β estimation per se. Rather, we use β to explore how the Q index relates to the statistical characteristics of the envelope.

3 Discussion and Conclusions

We have shown that foreshock waveforms often exhibit a distinct sawtooth-shaped envelope, characterized by multiple energy bursts and a prolonged decay. This pattern differentiates them from isolated events of similar magnitude and suggests a fundamentally different underlying process. A physical interpretation, consistent with the pre-slip model, is illustrated in Figure 4.

Prior to a large rupture, the impending rupture zone can be envisioned as a locked segment composed of multiple rate-weakening fault patches embedded within rate-strengthening surroundings (Figure 4a). Due to frictional heterogeneity, some patches approach failure earlier than others and may rupture. Two outcomes may follow. If adjacent regions are not critically stressed, the rupture remains confined, producing an isolated foreshock. In this case, the envelope function μ(t) decays after the main peak, reflecting local relaxation without broader fault activation. Alternatively, if the stress perturbation engages neighboring patches that are also near failure, a cascade may ensue, culminating in a large rupture that propagates across the fault system (Figure 4c). In this scenario, the resulting waveform exhibits a characteristic sawtooth-like envelope, reflecting successive bursts of seismic energy from spatially distributed failure. If the rupture continues to propagate and activates a large portion of the fault, it evolves into the mainshock of the sequence (Figure 4d). In this case, the envelope is dominated by high-amplitude fluctuations resulting from intense aftershock activity. The dense sequence of secondary events rapidly fills the signal, masking any sawtooth-like features that might have preceded the rupture. Thus, while moderate foreshocks may exhibit a distinctive waveform fingerprint, these signatures become increasingly difficult to resolve as the system transitions toward full rupture and seismic saturation.

Our results suggest that as the fault approaches failure, the envelope function μ(t) evolves into a sequence of discrete bursts separated by intervals where μ(t) ≃ μ_B. This pattern reflects the activation of isolated unstable patches within a locked region. The closer the fault is to global failure, the greater the number of critically stressed patches, and the more pronounced the sawtooth structure becomes. Following the mainshock, the redistribution of stress alters the spatial pattern of seismicity. Aftershocks tend to occur along the boundaries of the previously locked region, producing a sustained high plateau in μ(t) (Figure 4d). Based on these dynamics, we identify three characteristic envelope behaviors: (a) A single sharp peak in μ(t) followed by a rapid return to background levels (μ ≃ μ_B) indicates localized slip within a locked fault that remains far from global failure. (b) A sustained high plateau (μ_L ≫ μ_B) is consistent with intense aftershock activity following a mainshock, suggesting that no larger event is imminent in the near future. (c) A sawtooth-shaped envelope with multiple energy bursts separated by quiet intervals reflects a foreshock occurring within a fault region nearing global failure, where several unstable patches are interacting.

This framework accounts for the observed envelope of the 2011 M_w7.3 Tohoku foreshock, as well as the 2014 M_w6.2 Iquique foreshock, both of which show pronounced sawtooth patterns (black curves in Figures 2d–2f). These events likely ruptured within the interior of the locked region, where the interaction of multiple unstable patches generated distinct bursts of seismic energy. In contrast, foreshocks that occur near the boundary of the locked region show a different behavior. In these cases, stress is partly transferred outward, triggering normal aftershock sequences. This scenario is exemplified by the 2014 M_w6.7 Iquique foreshock (Figure 2f), where μ(t) decays to a low but nonzero plateau, indicating modest aftershock activity. At the same time, large isolated peaks reflect additional seismic bursts likely associated with foreshock dynamics.

A similar mixed pattern is observed in the 2019 Ridgecrest sequence. The M_w6.4 foreshock, located on a secondary fault adjacent to the eventual rupture zone, generated both a weak aftershock sequence and large-amplitude bursts linked to foreshock activity on the main fault. This results in a low but nonzero μ_L − μ_B value, superimposed with sporadic peaks (green curve, Figure 2h). Notably, a M_w5.4 foreshock occurred just 16 hr before the mainshock, within 2 km of its epicenter. Its envelope strongly resembles that of the M_w7.3 Tohoku foreshock, suggesting a common underlying mechanism (black curve, Figure 2h).

Our findings demonstrate that foreshock sequences can be distinguished from typical aftershock sequences by analyzing ground velocity envelopes. Foreshocks preceding large earthquakes exhibit anomalous sawtooth-shaped waveforms, unlike the plateau-like patterns observed in aftershocks. These oscillations between high and low envelope values indicate the presence of multiple rate-weakening patches approaching failure.

We introduced the quantitative index Q, based on waveform anomalies, to assess the likelihood that an earthquake is a foreshock. Applied to a global data set of M6+ events, the method successfully identified 10 foreshocks that preceded larger earthquakes within 10 days. The Q index is advantageous due to its straightforward formulation in terms of μ_B, μ_L, and n_obs, which can be estimated within minutes of the main envelope peak (see Text S2 in Supporting Information S1). This approach addresses limitations in prior methods, such as the one proposed by Gulia and Wiemer (2019), which relies on earthquake-specific adjustments to define a reference b-value, potentially introducing bias (E. Brodsky, 2019; Marzocchi et al., 2019). In contrast, the Q index is universally defined and largely unaffected by aftershock detection biases (de Arcangelis et al., 2018).

While the method is robust and operationally simple, several factors that may influence Q remain to be systematically explored. These include sensitivity to the frequency band used to compute the envelope, variations in station geometry and epicentral distance, source depth, and the finite extent of large ruptures. Although our results show consistency across stations and parameter choices, a more detailed analysis of these effects would strengthen the generality of the approach. We hope this work provides a foundation for future efforts to better constrain the physical and observational limits of waveform-based foreshock identification.

To advance short-term seismic hazard assessment, waveform-based methods like the Q index should be integrated with independent indicators of fault stress evolution, such as spatial and temporal variations in the b-value (Gulia & Wiemer, 2019). Large foreshocks often induce Coulomb stress changes that correlate with b-value anomalies (Nanjo, 2020), offering insight into the evolving state of the fault. For example, a sustained decrease in b-value was observed over several years prior to the 2011 Tohoku earthquake (Nanjo et al., 2012; Tormann et al., 2015), whereas an increase was documented near the mainshock hypocenter in the 2016 Kumamoto sequence (Nanjo & Yoshida, 2017). While the b-value captures longer-term stress accumulation, the Q index provides a near-instantaneous measure of waveform anomalies linked to cascading failure. Combined, these approaches offer complementary perspectives on fault loading and rupture potential. We suggest that future operational strategies for earthquake forecasting should explore multi-parameter frameworks that integrate waveform features with seismicity-based stress indicators.

Conflict of Interest

The authors declare no conflicts of interest relevant to this study.