# Do large earthquakes start with a precursory phase of slow slip?

Quentin Bletery, Jean-Mathieu Nocquet

► To cite this version:

Quentin Bletery, Jean-Mathieu Nocquet. Do large earthquakes start with a precursory phase of slow slip?. Seismica, 2025, 3 (2), 10.26443/seismica.v3i2.1383. hal-04886318

HAL Id: hal-04886318

https://hal.science/hal-04886318v1

Submitted on 16 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Do large earthquakes start with a precursory phase of slow slip?

Quentin Bletery ID * 1, Jean-Mathieu Nocquet ID 1,2

¹Observatoire de la Côte d'Azur, Université Côte d'Azur, IRD, CNRS, Géoazur, France, ²Université Paris Cité, Institut de Physique du Globe de Paris, CNRS, France

Author contributions: Conceptualization: Q.B., J-M.N.. Software: Q.B., J-M.N.. Formal Analysis: Q.B., J-M.N.. Writing - Original draft: Q.B.. Writing - Review & Editing: J-M.N.. Funding acquisition: Q.B., J-M.N..

**Abstract** In a recent publication, we showed that a stack of all GPS time series recorded before Magnitude ≥ 7.0 earthquakes suggests that large earthquakes start with a precursory phase of accelerating slow slip (Bletery and Nocquet, 2023). While no peer-reviewed comment or publication has formally contradicted this result, informal discussion has emerged on various platforms. We present here the different elements of discussion and address them through a series of tests. In particular, it has been proposed that correcting GPS time series from network common-mode noise makes the signal vanish. We confirm this result, but we show that this common-mode filtering procedure may inadvertently remove an existing tectonic signal. Moreover, the analysis of past records indicates that the likelihood that common-mode noise produces the signal we observe is well below 1 %. Additionally, we find that the signal is maximum at the location of impending earthquakes, and for a slip direction (rake angle) close to the one of upcoming events. The collective outcomes of these tests make very unlikely that the signal solely arises from noise. Even though the results of our tests do not irrefutably demonstrate the existence of a precursory phase of slow slip, they do support its existence. We hope that this study will motivate further work to provide a definite answer to the question of the tectonic origin of the observed signal and confirm or refute that large earthquakes start with a precursory phase of slow slip.

Production Editor:

Gareth Funning

Handling Editor:

Wenbin Xu

Copy & Layout Editor:

Théa Ragon

Signed reviewer(s):

Jeffrey Freymueller

Received:

June 5, 2024

Accepted:

December 11, 2024

Published:

December 31, 2024

# 1 Introduction

The search for precursory signals to large earthquakes has been a long-standing pursuit and the existence of such signals has generated a multi-decadal debate in the earthquake science community (Scholz et al., 1973; Geller, 1997; Kagan, 1997). A phase of precursory slip acceleration leading to the rupture is systematically seen in laboratory experiments (Ohnaka and Shen, 1999; Latour et al., 2013; Passelègue et al., 2017; Hulbert et al., 2019) and in dynamic models (Dieterich and Kilgore, 1996; Rubin and Ampuero, 2005; Kaneko et al., 2016). The duration of the precursory phases observed in experiments and models is on the order of a microsecond, but could become arbitrarily longer by considering heterogeneous faults (Lebihain et al., 2021), which are much more realistic than the homogeneous ones typically considered in experiments and models. With the development of geodetic and seismic instrumentation worldwide, observations of potential pre-seismic activity have been made on natural faults, suggesting the existence of a potentially observable precursory phase of slow slip on the fault preceding the rupture (Bouchon et al., 2011; Kato et al., 2012; Bouchon et al., 2013; Brodsky and Lay, 2014; Mavrommatis et al., 2014; Ruiz et al., 2014; Schurr et al., 2014; Bouchon et al., 2016; Radiguet et al., 2016; Ruiz et al., 2017; Soc-

quet et al., 2017; Ellsworth and Bulut, 2018; Tape et al., 2018; Bedford et al., 2020; Caballero et al., 2021; Beaucé et al., 2023; Martínez-Garzón and Poli, 2024). Nevertheless, those observations do not appear to be systematic and their causal relationship with the subsequent seismic events is not clear since the observations generally do not directly precede the earthquakes and similar ones are routinely made at times not preceding earthquakes (Schwartz and Rokosky, 2007; Gomberg et al., 2010; Obara and Kato, 2016; Bletery and Nocquet, 2020; Wallace, 2020; Behr and Bürgmann, 2021).

As observing a potential precursory slow slip acceleration seems out of reach at the scale of individual events, we conducted a global analysis of GPS displacement time series recorded before all large earthquakes (Bletery and Nocquet, 2023). For that purpose, we used all the available high-rate (sampled at 5 min) GPS time series, provided by the Nevada Geodetic Laboratory (NGL) (Blewitt et al., 2018), recorded in the 48 hours before Moment Magnitude ($M_w$) ≥ 7 earthquakes within a 500 km radius from the epicenter of the upcoming events (only excluding time series containing gaps or obvious offsets). We considered the hypocentral locations and focal mechanisms provided by SCARDEC (Valée and Douet, 2016) to compute the displacements $\vec{g}_{i,j}$ expected from a hypothetical precursory slow slip before each event $i$, at each GPS site $j$. We then calculated the dot product of each displacement measure-

<footnote>*Corresponding author: bletery@geoazur.unice.fr</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1** Stack of the dot product between the displacements expected from hypothetical precursory slip and the displacements recorded by GPS in the 48 hours preceding 90 Mw ≥ 7 earthquakes at 3,026 GPS stations (Bletery and Nocquet, 2023).</figcaption></figure>

ment $\vec{u}_{i,j}(t)$ (at station j, at time t before earthquake i) with the corresponding expected displacement $\vec{g}_{i,j}$, at each 5-min time step, and stacked all the obtained time series, resulting in a global stack S of 3,026 time series recorded before 90 earthquakes,

$$S(t) = \sum_{i=1}^{N_{eq}} \sum_{j=1}^{N_{st}(i)} \frac{\vec{u}_{i,j}(t) \cdot \vec{g}_{i,j}}{\sigma_{i,j}^2}, \quad (1)$$

where $\sigma_{i,j}$ is an estimate of the noise amplitude at each station (calculated as the L2-norm of $\vec{u}_{i,j}(t)$ from 48 to 24 hours prior to the events), $N_{st}(i)$ is the number of stations for earthquake i and $N_{eq} = 90$ is the number of earthquakes. The result of this dot product stack showed a subtle increase in the ~2 hours directly preceding the events (Figure 1), indicating a growing consistency between the recorded and the expected displacements as the faults approach failure, which we interpreted as indicative of accelerating precursory slip (Bletery and Nocquet, 2023).

Though this result has not been formally contradicted in a peer-reviewed comment or publication, several questions have emerged on informal platforms (e.g., Bradley and Hubbard, 2023a,b; Bürgmann, 2023; Voosen, 2023).

(1) How much does the uneven relative weight of the different events bias the stack?

(2) Does the observed signal arise from network-scale correlated (hereafter referred to as common-mode) noise?

(3) Does the observed signal originate from co-seismic contamination of the pre-earthquake GPS time series?

(4) May the signal be explained by foreshocks preceding some events?

In this study, we address these different questions through a series of tests with the objective of bringing new elements to the discussion on the origin of the observed signal (hereafter referred to as the signal).

## 2 Relative weights of the different stations and earthquakes in the stack

Given that the GPS stations are located at varying distances from the source of the impending earthquakes, the dot product between the GPS time series and the Green's functions have very different amplitudes (because the Green's functions have very different amplitudes). Moreover, the number of available observations drastically differs from one event to another. Consequently, the different events have different weights in the stack. One way to quantify the relative weights is to calculate the sum of the amplitudes of the Green's functions for the different events (normalized by the total sum),

$$\sigma_g(i) = \frac{\sum_{j=1}^{N_{st}(i)} |\vec{g}_{i,j}|}{\sum_{i=1}^{N_{eq}} \sum_{j=1}^{N_{st}(i)} |\vec{g}_{i,j}|} \quad (2)$$

$\sigma_g(i)$ provides an estimate of the intrinsic weight of earthquake i which is independent from the observa-

<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 2** a) Relative weights $\sigma_g(i)$ of the different earthquakes in the stack given by the sum of the amplitudes of the Green's functions (equation 2). b) Relative weights $\sigma'_g(i)$ accounting for the noise level (equation 3).</figcaption></figure>
tions. An alternative weight formulation is

$$ \sigma'_{g}(i) = \frac{\sum_{j=1}^{N_{st}(i)} |\vec{g}_{i,j}| / \sigma_{i,j}^{2}}{\sum_{i=1}^{N_{eq}} \sum_{j=1}^{N_{st}(i)} |\vec{g}_{i,j}| / \sigma_{i,j}^{2}}. \quad (3) $$

$\sigma'_{g}(i)$ is not independent from the observations but is closer to the weight in S as each time series is divided by the square of its estimated noise level $\sigma_{i,j}^2$ in the optimal stack (equation 1). We calculate the relative weight of each event in the stack using these two formulations (Figure 2).

The event with the largest weight ($\sigma_g = 18.9\%$ of the total, $\sigma'_g = 19.8\%$) is the 2010 M$_w$ 7.1 Baja California earthquake (Figure 2). The 2011 M$_w$ 9.0 Tohoku-Oki earthquake only arrives third with a weight $\sigma_g$ of 7.3% ($\sigma'_g = 9.4\%$). Therefore, the stack is not overly dominated by the Tohoku-Oki earthquake – nor any other earthquake – and the general shape of the stack is preserved when removing any individual event. Nonetheless, removing the 3 earthquakes with the largest weights – i.e. Baja California (2010, M$_w$ 7.1), Northern California (2014, M$_w$ 7.0), and Tohoku (2011, M$_w$ 9.0) – the signal strongly weakens. Note that this operation removes 837 observations (28% of the total), 39.3% of the expected signal (cumulative $\sigma'_g$) and the 3 best-recorded events (for which a precursory signal is most likely to be seen if it exists).

<figure><img src="image_5.png" /><figcaption>Figure 3 Cumulative amplitude of the Green's functions ($\sigma_g$) as a function of distance (in percentage of the total sum). The red dot indicates the weight of the data located closer than 200 km away from the epicenters relative to the total sum of the data located in a 500 km radius.</figcaption></figure>

One may think that only the nearest-field stations may contain information on a potential tectonic signal, as in the far field the amplitude of the Green's functions decreases with the square of the distance to the source (Mansinha and Smylie, 1971). This thinking neglects that the number of available stations also increases with the square of the distance to the source. As a result, when looking at the cumulative weight ($\sigma_g$) of the different observations as a function of their distance to the epicenter of the impending earthquakes (all events combined), we see that observations located more than 200 km away have a cumulative weight on the order of 30% of all observations within a 500 km radius (Figure

3). Note that the trend of the curve in Figure 3 suggests that even at distances larger than 500 km the cumulative weight of (noiseless) far-field observations will still increase more or less linearly with distance (as a consequence of the number of stations increasing at the same rate than the Green's functions amplitudes decrease). This means that, even though the signal-to-noise ratio of a potential tectonic signal strongly decreases with distance, making a potential signal invisible on individual far-field stations, this reasoning is not necessarily true when considering a stack of distant stations, and one should not assume that no tectonic signal can be visible in a stack of observations recorded farther than 200 km away from a potential source.

# 3 Sensitivity of the observed signal to common-mode filtering

## 3.1 Common-mode noise in GPS data

GPS time series contain noise correlated both in space and time, which overlaps with potential tectonic signals of interest. Over the past three decades, various approaches have been proposed to isolate and remove the regional common mode contribution in daily time series, where a single position is estimated from 24-hours long sessions (e.g., Wdowinski et al., 1997; Dong et al., 2006; Tian and Shen, 2016; Kreemer and Blewitt, 2021). Other research has focused on assessing time-correlated noise in individual residual time series (Zhang et al., 1997; Mao et al., 1999; Williams et al., 2004), while a few studies have examined both spatial and temporal noise characteristics simultaneously (e.g., Amiri-Simkooei, 2009; Gobron et al., 2024). The exact origin of both space and time correlated noise still remains unclear. Part of it arises from true physical motion of the ground induced by the response of the solid Earth to mass redistribution in continental hydrology, atmosphere and the ocean (e.g., Dong et al., 2002; Channard et al., 2018). In addition, mismodeling of the orbits or of the tropospheric delay has also been proposed to induce spatially correlated noise (e.g., Gobron et al., 2024). To our knowledge, there hasn't been similar systematic analysis of noise patterns in high-rate GPS time series, where positions are estimated at each measurement epoch. Unlike daily static analysis, which benefits from the averaging effect of numerous observations to estimate a limited number of parameters, high-rate GPS time series are directly affected by measurement noise, orbit mismodeling, atmospheric propagation delays, phase center variation correction errors, and multipath effects near the antenna. GPS satellites appear in the same part of the sky every sidereal day (23 hours, 56 minutes, 4 seconds). Since some GPS analysis errors are related to the receiver-satellite vector position, apparent displacement patterns in individual time series tend to repeat each sidereal day. This repetition property has been extensively used to post-process individual time series by removing these repeating patterns (e.g., Choi et al., 2004; Larson et al., 2007). Although this approach has proven to be efficient for periods of tens of seconds to tens of minutes, it is less clear how efficient
<figure><img src="image_6.png" /><figcaption>**Figure 4** a) Dot product stack of the raw time series for the 31 earthquakes recorded by at least 10 stations located farther than 200 km away from the epicenters. b) Dot product stack applied to the time series obtained after removing translational common modes for the same global stack. c) Dot product stack of the translational common modes alone.</figcaption></figure>

it is to remove hour-to-day long components of errors and whether it reduces the spatially correlated component of noise.

In our original study, we did not evaluate common mode errors before applying the multi-earthquake stacking procedure. Our strategy was to extract potential tectonic signal – aligning with displacements expected from hypothetical precursory slip – from the raw data. This strategy was intended to minimize subjective post-processing choices potentially biasing the analysis. Using raw time series allowed us to include all observations and treat each time series uniformly. Since there is no reason that common-mode noise aligns with displacements expected from fault slip, we assumed that common-mode noise will cancel out when stacking numerous earthquakes that have been recorded on distant local networks at distant times. The fact that the different earthquakes have different weights in the stack (Figure 2) makes the aforementioned assumption potentially questionable (Bradley and Hubbard, 2023a,b). We investigate, here, the potential role of common-

mode noise in the signal we observe in Figure 1.

## 3.2 Removing translational common modes before stacking makes the signal vanish

One way to evaluate network-scale correlated noise is to calculate a translational common mode as the mean of the GPS time series which presumably do not contain signals of interests (i.e. time series recorded at locations distant from the hypothetical sources). Given the heterogeneity of the datasets available for the different events, this is not possible for all the earthquakes. Restricting the dataset to the 31 events for which at least 10 stations are available more than 200 km away from the epicenter, we first verify that the stack presents a similar pattern to the one including all the earthquakes (Figure 4.a compared to Figure 1). We then evaluate a translational common mode, calculated as the mean of the time series recorded more than 200 km away from the epicenter of these 31 events. We remove this common mode from the time series and calculate
<figure><img src="image_7.png" /><figcaption>**Figure 5** a) Dot product stack of noiseless synthetic time series. b) Dot product stack of noise time series.</figcaption></figure>

the stack again. We verify that, as discussed on informal platforms (Bradley and Hubbard, 2023a), after such common-mode filtering, the signal in Figure 1 and 4.a can no longer be seen (Figure 4.b). Moreover, as also discussed on informal platforms (Bradley and Hubbard, 2023b), prescribing the per-earthquake estimated common modes as pseudo-observations to all time series, we obtain a stack time series in which a similar signal appears (Figure 4.c).

Figure 4 suggests that the origin of the signal is not tectonic but rather originates from network-scale correlated noise. Two questions arise at this stage. How likely is network-scale correlated noise to produce such a signal right before the earthquakes? And, is it possible that removing a translational common mode inadvertently removes tectonic signal? To address the first question we will quantify such likelihood in section 4. To address the second question, we can impose a synthetic precursory signal mimicking the proposed one, add noise, and test how the common-mode filtering procedure performs at separating the imposed tectonic signal from noise.

## 3.3 Is removing a translational common mode an efficient way to separate noise from weak tectonic signal?

We impose a growing slip – following the exponential fit of the global stack in Figure 2 of Bletery and Nocquet (2023) – on $1 \times 1$ km faults centered on the hypocenter of the 31 events for which evaluating a common mode is feasible. We calculate the synthetic displacements at the different GPS sites corresponding to the imposed slip and add noise. Because GPS satellites are seen at the

same location in the sky every sidereal day (i.e., every 23 h 56 min 4 s) (see section 3.1), spatial and temporal patterns of noise in high-rate GPS time series tend to repeat from one sidereal day to the next. In order to mimic as closely as possible the network-scale correlated structure of the noise in the day preceding the earthquakes, we therefore use the time series recorded from 48 to 24 h (minus 1 sample) before the earthquakes as a realistic noise including realistic network-scale correlations. We then calculate the stack before and after removing translational common modes, as defined in the previous paragraph. Since we imposed the tectonic signal and the noise, we know the target signal that an ideal noise filtering procedure should find. To visualize this target, we separately calculate the stack of the noiseless synthetic time series (Figure 5.a) and the stack of the noise time series (Figure 5.b).

Comparing the stacks obtained before and after removing the common modes and the stack of the common modes alone to the target stack (Figure 5.a), we find that the common-mode filtering performs poorly at separating the imposed signal from noise (Figure 6). As expected (since we imposed it), an exponential-like signal appears in the stack before removing the common modes (Figure 6.a). More surprisingly, the signal can no longer be seen in the stack after removing the common modes (Figure 6.b). Even more surprisingly, an exponential-like signal appears in the stack of the common modes alone (Figure 6.c) in a similar way as in Figure 4.

In the last 2 hours before the events, the misfit of the imposed signal with the denoised stack (Figure 6.b) is 20 % larger than with the stack of the common modes alone (Figure 6.c). This highlights that separating po-
<figure><img src="image_8.png" /><figcaption>**Figure 6** Same as Figure 4 for synthetic time series made from imposed tectonic signal plus network-scale correlated noise. The orange curve super-imposed on the 3 plots represents the stack obtained from the noiseless imposed signal (Figure 5.a), i.e. the target of the denoised stack.</figcaption></figure>

tential weak tectonic signal from network-scale correlated noise is a complex problem and that the basic translational common-mode filtering procedure – consisting in removing the average of time series recorded more than 200 km away from the epicenters – may improperly remove tectonic signal.

On the other hand, the results of the presented test are highly sensitive to the considered noise. Using different time windows as noise gives different pictures. In the general case, removing translational common modes makes the imposed signal more visible, but the improvement is not systematic. This highlights that the space-time structure of the noise is complex and that it is challenging to filter it without altering a potential weak tectonic signal.

## 3.4 Structure of the correlated noise in the 48 hours preceding the events

In order to explore the evolution of the structure of the network-scale correlated noise in the 48 hours preced-

ing the events, we calculate the cross-station dot product for each pair of stations before each earthquake at each increment of time and represent them as a function of distance (by bins of 10 km) and time windows of 2 hours (Figure 7). The cross-station dot product can be seen as a measure of the correlated noise. It is expected to be larger when stations are close to each other (left part of the curves) and then to stabilize in a plateau, as can be seen in the average on the [-48, -2] h time window preceding all the events (blue curve in Figure 7).

We can infer from Figure 7 that the network-scale correlated noise – i.e. the common-mode noise – is larger in the 2 hours preceding the earthquakes (red curve) than the average of the 2 days before (blue curve). We also see that it is significantly larger than the amplitude of the correlated noise 24 hours before (orange curve), suggesting that the main source of this elevated correlated noise is likely not orbital modeling errors (which tend to repeat from one day to the next). Moreover, other 2-h-long time windows exhibit larger cross-station correlations in the 48 hours preceding the events
<figure><img src="image_9.png" /><figcaption>**Figure 7** a) Cross-station dot-product as a function of distance for different 2 h time windows ranging from 48 to 46 h before the events for the first one, to 2 to 0 h for the last one (grey curves). The [-2, 0] h time window is highlighted in red, the [-26, -24] h one in orange, and the average for the [-48, -2] h time window in blue.</figcaption></figure>

than the [-2, 0] h one. This means that the signal in Figure 1 does not correspond to an exceptional pattern in the structure of the correlated noise as inferred from the cross-station dot product.

As it is difficult to isolate potential tectonic signal from network-scale correlated noise and to conclude on the origin of the signal observed in the last hours of the stack based on the study of the structure of the correlated noise, we will focus, in the next section, on the statistical significance of the signal.

# 4 Statistical significance of the signal

How likely is network-scale correlated noise to produce the signal we observe in Figure 1? This likelihood may be assessed by estimating how frequently similar patterns emerge from noise and if they would emerge assuming different source locations or different focal mechanisms.

## 4.1 The signal points to the time of the upcoming earthquakes

In Bletery and Nocquet (2023), we provided a first estimation of the likelihood that the signal originates from network-scale correlated noise by calculating the stack for a large number of time series recorded at random times on the stations considered in our original analysis. In each case, we calculated the ratio *r* between the last point of the stack moving average (with a moving window of 1 h 50 min) and its maximum in the 46 preceding hours (Figure 8.a). We found a value of *r* equal or larger to the one obtained using the time series preceding the earthquakes (*r* = 1.82) in 0.3% of the cases (Figure 8.b). We also counted the number *n* of monotonically increasing points at the end of the stack time

series and found values equal or larger to the one using the time series preceding the earthquakes (*n* = 23) in 0.8% of the cases (Figure 8.c). Combining the two, we found that *r* ≥ 1.82 and *n* ≥ 23 for 0.03% of the drawn time windows, providing a rough estimate of the likelihood that such a signal arises by chance from noise (see supplementary material of Bletery and Nocquet (2023) for details).

To further evaluate the probability that correlated noise in individual event stacks constructively sum up to produce a signal similar to the one we observed, we perform a complementary test by simulating 100,000 surrogates of stack times series for each earthquake. For that purpose, we randomly shuffle the phase of individual earthquake stack time series, preserving their Fourier amplitude. This allows us to synthetically simulate 100,000 stack time series for each earthquake that share the same characteristics than the original ones. We then calculate the 100,000 associated global stacks. We find values very consistent with the previous test: *r* ≥ 1.82 in 0.2% of the cases, *n* ≥ 23 in 0.9% of the cases, and the 2 combined in 0.02% of the simulated samples (Figure 8.d-e).

These two tests consistently indicate that network-scale correlated noise may coincidentally sum up constructively to produce a signal similar to what we observe but the likelihood of such a thing to happen precisely at the time we observe it is on the order of 0.03 % (0.3 % if we only consider *r*). We emphasize that these two tests provide statistics that take into account the uneven relative weight of the different events and network-scale correlated noise, overall indicating very low likelihood that the signal randomly arises from (common-mode) noise.
<figure><img src="image_10.png" /><figcaption>**Figure 8** a) Sketch illustrating the calculation of the ratio *r* between the last point of the moving average and the maximum of the moving average in the [-48,-2] h time period. b) Histogram of the ratio *r* between the last point of the moving average and its maximum in the preceding 46 hours for 100,000 random noise time windows (Bletery and Nocquet, 2023). c) Histogram of the number *n* of monotonically increasing points at the end of the moving average for the same 100,000 random noise time windows (Bletery and Nocquet, 2023). d) Histogram of *r* for 100,000 random surrogates of dot product stacks. e) Histogram of *n* for the same 100,000 random surrogates of dot product stacks. The vertical red lines show the values obtained for the moving average preceding the earthquakes.</figcaption></figure>

## 4.2 The signal points to the location of the upcoming earthquakes

The likelihood that the signal originates from network-scale correlated noise may also be assessed by its sensitivity to the spatial structure of the recorded displacements: if network-scale common noise dominates the recorded displacement time series, then randomly per-

muting the Green's functions (among GPS sites that recorded the same earthquake) should yield similar stacks. We test this idea and randomly shuffle the Green's functions associated with the different time series, earthquake by earthquake. We then stack together the stacks obtained for the different earthquakes and calculate the ratio *r* as previously defined. On 100,000 random permutations of the Green's functions, we find
<figure><img src="image_11.png" /><figcaption>**Figure 9** a) Histogram of *r* for 100,000 random permutations of the Green's functions. b) Same as a) excluding the Tohoku-Oki earthquake. The orange and red lines respectively show the median of the distributions and the value with the correct Green's functions.</figcaption></figure>

a median value of *r* of 1.23 (Figure 9.a). This value is very high and suggests that a significant part of the signal may be related to common-mode noise. Nevertheless, we find that *r* ≥ 1.79 (value with the correct Green's functions excluding events recorded by only 1 station, for which shuffling the Green's functions is not possible) for only 6.5 % of the permutations. This last number may be seen as an alternative estimate of the likelihood that the signal arises solely from common-mode noise based on the spatial structure of the signal.

<figure><img src="image_12.png" /><figcaption>**Figure 10** Ratio *r* between the last point of the moving average and its maximum on the [-48,-2] h time window for stacks calculated considering sources on a grid of 81 locations spaced by 50 km and centered on the correct source locations. We find that *r* is maximum at the correct source locations (denoted by the red star).</figcaption></figure>

Since the Tohoku-Oki earthquake was preceded by a significant foreshock (51 hours before the mainshock), a possible afterslip signal following this foreshock may arguably bias the Tohoku stack. We therefore reproduce the test above excluding the Tohoku event. This changes the value of *r* with the correct Green's functions to 1.42

<figure><img src="image_13.png" /><figcaption>b) Distribution of r obtained by shuffling the Green's functions excluding Tohoku</figcaption></figure>

and the median shuffling the Green's functions to 0.88 (Figure 9.b). Overall, we find that *r* ≥ 1.42 for 8 % of the random permutations, confirming that it is unlikely that the spatial structure of the signal emerges solely from noise.

Another way to assess whether the signal is most likely resulting from network-scale correlated noise or from a tectonic process related to the upcoming earthquakes is to alter the location of the sources before calculating the Green's functions. Moving the sources 100 km away in the east, west, north and south directions, the obtained stacks do not show a signal similar to the stack calculated considering the correct locations (Figure 11). We generalize the test and calculate *r* considering locations on a 400 × 400 km grid spaced by 50 km. We find that *r* is maximum at the actual location of the earthquakes (Figure 10), strengthening the idea that the signal originates from precursory slip in the direct vicinity of the hypocenters of the impending earthquakes.

## 4.3 The signal points to the slip direction of the upcoming earthquakes

Another way to assess how likely the signal we observe in Figure 1 is to be related to the upcoming earthquakes is to perturb the focal mechanism of the earthquakes in the calculation of the Green's functions before computing the stacks. When perturbing the rake angle $\lambda$ by large values ($\Delta\lambda \in [-180^\circ, -90^\circ, 90^\circ]$) the signal completely vanishes (Figure 12.a-c). We generalize the test by calculating *r* for rake perturbation increments of 10° from -180° to 170°. We find that *r* is large ($r \ge 1.5$) only for small perturbations of the rake angle ($|\Delta\lambda| \le 30^\circ$) – and *r* < 1.5 for $|\Delta\lambda| > 30^\circ$, further suggesting that the signal is related to the upcoming earthquakes.
<figure><img src="image_14.png" /></figure>

**Figure 11** Stacks obtained using Green's functions calculated considering sources moved 100 km away from the correct source locations in the east (a), west (b), north (c) and south (d) directions to be compared to the stack considering the correct locations (e). The red curves show the moving average of the different stacks. The black dashed lines represent the 0 line and the maximum of the moving average on the [-48,-2] h time window.
<figure><img src="image_15.png" /><figcaption>**Figure 12** a) Stack obtained after perturbing the rake of the earthquakes by $\Delta\lambda = -180^\circ$ (or $180^\circ$) from the catalog value. b) Same as (a) for $\Delta\lambda = -90^\circ$. c) Same as (a) for $\Delta\lambda = 90^\circ$. d) Value of $r$ obtained for perturbations of the rake ($\Delta\lambda$) going from $-180^\circ$ to $170^\circ$.</figcaption></figure>

The collective outcomes of the 5 tests presented in this section outline that, though subtle and not robust to translational common-mode filtering, the signal points to the time, the location and the mechanism of the impending earthquakes with a high statistical significance. This makes the tectonic origin of the signal more likely than network-scale correlated noise which has no reason to point to the time, location and mechanism of the events.

# 5 Contamination by co-seismic signal?

It has also been suggested that the observed signal could be an artifact resulting from the strategy used in the GPS analysis, which would tend to bias pre-earthquake positions by a fraction of the subsequent co-seismic offsets. The time series we used were processed, and gra-ciously made available to the community, by NGL (Ble-witt et al., 2018). They result from Precise Point Po-
<figure><img src="image_16.png" /><figcaption>**Figure 13** Dot-product stack replacing the Green's functions by the co-seismic offsets (excluding the 20 stations for which co-seismic offsets are unavailable).</figcaption></figure>

sitioning (PPP) kinematic analysis using GipsyX (http://geodesy.unr.edu/gps/ngl.acn.txt). The position is determined using carrier phase measurements decimated to every 5 minutes. The 5-minute pseudorange is computed by averaging the higher rate (typically 30 s) points against the carrier phase, but is only effectively used to enable carrier phase ambiguity resolution, after which the pseudorange contribution to the solution is completely negligible. Independently of whether an earthquake happened or not, the positions are formally correlated because of common parameters in the least-squares estimation. Common parameters are zenith troposphere, two tropospheric gradients and the carrier phase ambiguity. However, ambiguity resolution effectively breaks the covariance between positions and carrier phase ambiguity parameters because they become perfectly known, meaning that these correlations exist but are independent of the actual station motion and are independent of earthquakes (Geoffrey Blewitt, personal communication).

Station coordinates are estimated as random walk with a very large process noise, so that, effectively, there is no forced correlation between 5-minute epoch estimates, allowing station coordinates to “jump” to completely different values. The filter is first run forward in time, being blind to the future. Then, the filter takes the final estimated parameter state, and moves backward in time. In the NGL analysis, the process noise is set to $1 \text{ m.s}^{-1/2}$. This means that the a priori constraint controlling the change of position between adjacent epochs is ~17 m for 5-min samples. At any given epoch, the estimate from the next future epoch influences the current epoch with an a priori sigma of 17 meters when estimating the position using least-squares (Geoffrey Blewitt, personal communication). This constraint assigns a weight to the smoothing that is many orders of magnitude smaller than the weight of the carrier phase measurements which have precision of ~1 cm. Although such a loose smoothing constraint is likely too small to cause co-seismic offsets to bias positions before the earthquakes, the impact of this parameter could be tested in future studies.

Even though the NGL analysis is expected to prevent any co-seismic contamination, we investigate this possibility by replacing the Green's functions (computed considering a point source) originally used to calculate the

stack by the co-seismic offsets (the difference between the first measurement after the event and the last before the earthquake) recorded at each station (Figure 13). The idea behind this test is that if the signal in the last 2 hours is an artifact of co-seismic leakage, the artifact should be strongly correlated with the recorded co-seismic offsets (from which it presumably leaks from) and should appear stronger when taking the dot product with the co-seismic offsets than with the Green's functions.

We see from Figure 13 that the signal is not stronger when replacing the Green's functions by the co-seismic offsets in the global stack. Note that the reason the signal disappears is most likely because the recorded offsets are dominated by noise for many of the stations. Nevertheless, if the signal resulted from a problem of filtering leakage, even noise would leak, making the signal more apparent when replacing the Green's functions by co-seismic offsets consistently determined from the data set used in the pre-earthquake analysis. This indicates that the signal is not particularly correlated with the co-seismic offsets, suggesting that the signal is unlikely to result from co-seismic contamination. Nevertheless, the aforementioned test does not provide a definite answer to the question of a possible co-seismic contamination of the pre-seismic time series as one could imagine that more complex and subtle contamination processes would not necessary result in high correlations with the recorded co-seismic offsets. For instance, centimeter level co-seismic offsets at a few sites from the global tracking network could induce biases in orbit/satellite clock determination that would in turn leak into positions, possibly as a long wavelength common mode motion. In all the presented tests and in our original study, we rely on the only globally homogeneous GPS dataset made available by the Nevada Geodetic Laboratory. Independent GPS analyses would be informative to infer the sensitivity of potential pre-earthquake signals to different GPS analysis strategies, such as the possible impact of co-seismic static and dynamic motion of ground stations used to determine the satellite orbit and clock products.
<figure><img src="image_17.png" /><figcaption>**Figure 14** Same as Figure 7 for the Tohoku earthquake alone.</figcaption></figure>

## 6 The case of the Tohoku-Oki earthquake

In our original study, we treated the Tohoku-Oki earthquake as a special event. It was, by far, the largest event in the database ($M_w$ 9.0) and – even though it was not the one with the largest weight in the stack (Figure 2) – it was the event for which we had the largest number of observations ($N_{st}$ = 355). Therefore, the Tohoku-Oki earthquake was the one event for which we were hoping that observing a signal at the scale of an individual earthquake could be possible. When looking at the stack obtained in the 24 h preceding the Tohoku-Oki earthquake, we observed an unexpected seemingly-periodic signal possibly super-imposed on an exponential-like one (Bletery and Nocquet, 2023).

We quantified how exceptional the periodicity of this signal was by calculating the misfit reduction provided by a sinusoidal function defined as $y = A \sin(t + \phi) + B$. The obtained misfit reduction appeared to be exceptional compared to stacks calculated at other times and considering other source locations (Bletery and Nocquet, 2023). We realized that this exceptional misfit re-

duction was not due to an exceptional periodicity but to a large value of $B$, likely due to afterslip that developed between the 2011 March 9 foreshock and the main-shock. Estimating the misfit reduction arising from the periodic oscillation alone, the periodic signal observed before the Tohoku earthquake does not appear to be unique. This invalidates the interpretation we made of this seemingly-oscillatory behavior as a potential precursor signal and rather suggests that the oscillations originate from network-scale correlated noise.

This also raises the question of the origin of the signal we observed in the final hours before the Tohoku earthquake. To investigate this question we apply the cross-station dot product calculation (see section 3.4) to the data recorded before the Tohoku event alone. It reveals a different picture than in the global case (Figure 14). The cross-station dot product appears larger in the last 2 hours before the event than in any other time window in the 2 days before, including the one 24 hours before (Figure 14). This suggests large common-mode noise at that particular time, which we do not observe for any other event and which we do not observe – to this point – on average (Figure 7) despite the effect of the

<figure><img src="image_18.png" /><figcaption>**Figure 15** Same as Figure 13 for the Tohoku earthquake alone (orange) compared to the original stack for the Tohoku earthquake (blue) excluding the 2 stations for which co-seismic offsets are unavailable.</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>**Figure 16** a) Relative weights $\sigma_g(i)$ of the different earthquakes for the updated stack (equation 2). b) Relative weights $\sigma'_g(i)$ for the updated stack (equation 3). Light blue slices indicate events added in the update.</figcaption></figure>

<figure><img src="image_20.png" /></figure>

Tohoku data (included in Figure 7). This behavior could be indicative of (1) an unfortunate large common-mode noise (likely not due to orbital miss-modeling as the [-26,-24] h time window does not exhibit the same pattern), (2) co-seismic contamination of the pre-Tohoku time series, but (3) would also be consistent with our original interpretation of precursory slip on the fault area surrounding the hypocenter of the upcoming Tohoku earthquake.

To investigate the second possibility, we calculate the stack replacing the Green's functions by the co-seismic offsets as in the global case, for the Tohoku event alone. The result is more ambiguous than in the global case, with a stack obtained with the co-seismic offsets very similar to the original one but not exhibiting a stronger signal (Figure 15). This is somehow to be expected as, in this case, the co-seismic offsets are fairly similar to the Green's functions – given the magnitude of the event ($M_w = 9$), the co-seismic signal is many times larger than the noise at every station – and not particularly indicative of co-seismic contamination since the signal does not appear more clearly than in the original stack (Figure 15).

Overall, it is difficult to conclude whether the final positive increase before the Tohoku-Oki event is due to a precursory process, common-mode noise or co-seismic contamination. A recent study applying the stacking procedure we proposed to tilt records reports no evidence of slow slip preceding the Tohoku-Oki earthquake, indicating that if there was one, its cumulative moment magnitude was below 6.4 (Hirose et al., 2024). The level of noise (correlated at the scale of one network) makes it difficult to analyze the stacks obtained for individual events. Therefore, even though the presented tests suggest the possibility of the existence of precursory signal preceding the Tohoku earthquake and encourage further work in that direction, we do not conclude on the specific case of this event.

# 7 Update on recent earthquakes

We update the stack in Bletery and Nocquet (2023) by adding GPS time series recorded before recent earthquakes (Figure 16). The updated stack includes time series recorded on 5,015 stations before 109 earthquakes

(against 3,026 stations and 90 events in the original dataset). Among the added events, 4 have a significant weight: 2 earthquakes that happened offshore Honshu (Japan) in February 12 ($M_w$ 7.2) and March 20 ($M_w$ 7.1) 2021, the 2023 $M_w$ 8.0 Kahramanmaraş earthquake (Turkey) and the 2024 $M_w$ 7.6 Noto earthquake (Japan).

Given the proximity of station J253 to the Noto earthquake hypocenter, using hypocenter locations provided by different agencies leads to drastic changes in the direction of this station's Green's function and consequently – given the large amplitude of this Green's function – to significant changes in the global stack itself. Because of the sensitivity of the stack to location errors for this particular data point, we remove J253 from the stack. The shape of the updated stack exhibits large high-frequency fluctuations (such as the original one) but still highlights a positive increase at the end of the time series with a duration similar to the original stack (Figure 17.a). In fact, even though a high-frequency fluctuation makes the stack go down in the last minutes before 0, the *r* ratio increases to 2.1 (Figure 17.b) compared to 1.82 in Bletery and Nocquet (2023) (*r* = 2.06 if we do not remove J253). Using a time window of 3 hours gives *r* = 2.46 (Figure 17.c).

Even though the signal is arguably not as visually impressive as in Figure 1 because of a high-frequency negative trend in the minutes preceding the events, the positive trend in the previously identified time window (1 h 50 min) is actually strengthened by the addition of the recent events (*r* = 2.1). This result strongly encourages regular updates of the stack as newly-acquired data preceding large events become available. As exemplified by the 2024 Noto event – the new best-recorded event in terms of number of observations (695 stations) and of weight of the Green's functions (Figure 16) –, earthquakes to come will likely bring more and more information that will eventually confirm or refute the existence of an average slow slip acceleration leading up to large earthquakes.
<figure><img src="image_21.png" /><figcaption>**Figure 17** a) Updated stack including recent earthquakes. b) Moving average of the updated stack using the same time window as in Bletery and Nocquet (2023): 1 hour 50 minutes. c) Same as b) using a time window of 3 hours.</figcaption></figure>

# 8 Discussion

## 8.1 Responses to the questions asked by the community

In the introduction, we identified 4 questions that were several times asked by colleagues after the publication of our original study. We address them below.

### 8.1.1 How much does the uneven relative weight of the different events bias the stack?

The uneven relative weight of the different events is at the very basis of our stacking approach. The dot product with the Green's functions gives a natural weight to the observations that is suitable to extract weak signal from noise in an optimal stack. As illustrated in section 2, this results in some events counting significantly more than others in the stack. If the data were all independent from each other, this would not constitute a problem. However, since GPS time series are correlated in space and time at the scale of a regional network, this potentially gives a lot of weight to network-scale cor-

related noise recorded before events that have a large weight in the stack. A first indication that the signal we observe is not the result of a bias caused by the uneven relative weight of the different events is that adding recent events – some of which having a very large weight (Figure 16) – strengthens the significance of the signal (Figure 17).

### 8.1.2 Does the observed signal arise from network-scale correlated noise?

Can the signal we observe be due to an unfortunate combination of common-mode noise – aligned with the direction of the Green's functions – recorded before events that have a large weight in the stack? A first quick answer to this question is yes, as removing translational common modes – estimated as the mean displacement time series recorded by stations located more than 200 km away from the potential sources – removes the observed signal (Figure 4). Nevertheless, because the number of observations increases with distance at the same rate as the amplitude of a tectonic signal is expected to decrease, non-negligible tectonic signal con-
tribution in the stack may come from far-field stations (Figure 3). Consequently, the assumption behind the estimation of common-mode noise that far-field stations do not contain tectonic signal may be inaccurate. Consistently, we find that when imposing a synthetic signal, the aforementioned common-mode removal procedure inadequately identifies tectonic signal as noise – and noise as signal – (Figure 6), highlighting that there is a definite possibility that a real precursory signal would vanish after removing common modes estimated this way.

Moreover, we find – through 5 independent tests accounting for both the uneven relative weights of the events and common-mode noise – that, though subtle and not robust to common-mode filtering, the signal points to the time, location and slip direction of the upcoming events with a high statistical significance (section 4). This finding is a strong indication that the signal is unlikely to originate solely from network-scale correlated noise.

### 8.1.3 Does the observed signal originate from co-seismic contamination of GPS time series?

An alternative hypothesis that would explain the space-time structure of the signal (pointing to the time, location and mechanism of the events) would be that the signal originates from co-seismic contamination of the pre-earthquake data. A quick estimation of the potential bias in the GPS analysis of NGL points to a negligible effect. Nevertheless, a controlled experiment (manually moving an antenna) would be worth performing to rigorously estimate this bias. Moreover, we find that replacing the Green's functions by the co-seismic offsets (that the signal would presumably leak from) in the stack calculation does not strengthen the signal (and even makes it vanish), suggesting that the signal is not an artifact of co-seismic leakage. In all the presented tests, we rely on the only globally homogeneous GPS dataset made available by the Nevada Geodetic Laboratory. Independent GPS analyses would also be informative to infer the sensitivity of potential pre-earthquake signals to different GPS analysis strategies.

### 8.1.4 May the signal be explained by foreshocks preceding some events?

Earthquakes are known to occur in clusters (e.g., Helmstetter and Sornette, 2003). Consequently, large earthquakes are often preceded by foreshocks (e.g., Jones and Molnar, 1979; van den Ende and Ampuero, 2020; Moutote et al., 2021). Comments arising from the community suggested that the signal we observe could be due to such foreshocks (e.g., Voosen, 2023). In order to produce the signal we observe, the cumulative seismic moment of these events should correspond to an equivalent magnitude of 6.3. If, as Figure 9 suggests, part of the signal is due to common-mode noise, the cumulative moment could be reduced but could not go below an equivalent magnitude of 5.6. Foreshocks of such magnitude would clearly be seismically visible and catalogued as such, meaning that if they were at the origin of the signal, we should record, on average, a $M_w \ge 5.6$ seismic

event in the 2 hours preceding each $M_w \ge 7$ earthquake. Since this is clearly not the case, we do not believe that foreshocks are a plausible explanation for the signal we observe.

## 8.2 Additional questions

### 8.2.1 Have we used relevant statistical indicators?

The statistical tests we performed – both in Bletery and Nocquet (2023) and this study – mainly rely on two indicators: $r$ and $n$. Both of these indicators are calculated on a moving average using a moving window of 1 h 50 min. This time window is arbitrary and different ones would give different statistics. We see, for instance, that applied to the stack updated with the recent earthquakes, $n$ (the number of monotonically increasing points at the end of the moving average) is drastically reduced because of a high-frequency negative trend directly preceding the ruptures (Figure 17.b). Changing the moving window drastically changes $n$ (Figure 17.c). This illustrates that $n$ is probably not the most relevant statistical indicator. The ratio $r$ (that we used the most) between the last point of the moving average and its maximum on the rest of the time series is a lot more stable: changing the moving window does not change much $r$. The $r$ indicator is also a fairly intuitive proxy for a signal to noise ratio: the last point of the moving average is nothing more that the mean displacement in the last 1 h 50 min and the maximum of the moving average in the preceding 46 hours is a good measure of the noise fluctuations filtered at the period of interest. We believe $r$ is a reasonable statistical indicator, but it will be interesting to reproduce the statistics we obtained using other statistically-relevant indicators.

### 8.2.2 What is the effect of the point source approximation?

In Bletery and Nocquet (2023), we considered point-source-like sources in the calculation of the Green's functions (in practice, very small $1 \times 1$ km finite faults). The rationale behind this choice was that (1) models of earthquake nucleation usually involve a portion of the fault which is much smaller than the subsequent earthquake area, and (2) the point-source approximation allowed us not to have to make any a priori assumption on the extent of potential pre-slip faults. After a careful selection of the known or most probable nodal plans, we test the influence of considering extended sources of different lengths $L$ and widths $W$. We find that the result of the stack is fairly insensitive to the size of the considered source (Figure 18). Nevertheless, the ratio $r$ consistently increases with larger fault areas: $r = 1.89$ for ($L = 10$ km, $W = 10$ km), $r = 1.93$ for ($L = 20$ km, $W = 20$ km), $r = 1.98$ for ($L = 50$ km, $W = 20$ km), and $r = 2.02$ for ($L, W$) corresponding to the extent of the co-seismic rupture (following the scaling law empirically derived by Wells and Coppersmith, 1994). One may interpret this observation as suggesting that precursory slip occurs on large fault portions – possibly of size equivalent to the final rupture (see section 8.2.4) –
<figure><img src="image_22.png" /><figcaption>**Figure 18** Stack obtained considering extended sources of length *L* and width *W* in the calculation of the Green's functions: a) *L* = 10 km, *W* = 10 km, b) *L* = 20 km, *W* = 20 km, c) *L* = 50 km, *W* = 20 km, d) *L* and *W* corresponding to the co-seismic slip areas of each event based on the empirical scaling laws derived by Wells and Coppersmith (1994).</figcaption></figure>

but we believe the changes in the stack are too small to
support this interpretation.

8.2.3 Is precursory accelerating slow slip system-
atic?

Assuming the signal we observe is generated by an ac-
celerating slow slip, is this behavior systematic or is it
resulting from only a few events? Given that the signal
we observe is at the very limit of the detection threshold

in the global stack, we only have access to the average
behavior prior to all the events. Inferring precursory
signal at the scale of individual events – or even subsets
of events – is out of reach. Therefore, we cannot con-
clude on whether the proposed signal originates from
all events or a specific subsets of them.
### 8.2.4 Does precursory slip depend on magnitude?

A natural related question is whether or not the amplitude of the proposed precursory signal scales with magnitude, as laboratory experiments suggest (e.g., Acosta et al., 2019). It seems plausible that some kind of scaling exists – as it would seem illogical that a magnitude 1 event produces an accelerating slow slip of equivalent magnitude 6.3 – but here again, the available data does not allow us to answer the question.

## 8.3 Perspectives

The most important pending question is the possible influence of network-scale correlated noise in the signal we observe. The translational common-mode estimation presented in this study is only one among many existing approaches to mitigate noise in GPS time series. Alternative – more sophisticated – approaches such as Independent Component Analysis (ICA) or variational bayesian ICA (Gualandi et al., 2016) will be interesting to apply. Regular updates of the stack including events to come will also be informative and, provided enough time, will eventually confirm or refute the existence of the signal. Other perspectives include reproducing our results using independent GPS solutions (the only global one presently available is the NGL one), analyzing smaller magnitude events, and looking at other types of data. For instance, one would expect that a slow slip acceleration generates an increase in micro-seismic activity as is observed during weeks-long slow slip events (Schwartz and Rokosky, 2007; Gomberg et al., 2010; Obara and Kato, 2016; Bletery and Nocquet, 2020; Wallace, 2020; Behr and Bürgmann, 2021). Analyzing the evolution of micro-seismic noise recorded by seismic stations located in the vicinity of the source of large earthquakes in the hours preceding their initiation could reveal crucial complementary information on the nucleation phase of these events.

## 9 Conclusion

We built on the global analysis of GPS time series preceding large earthquakes that highlighted an average growing displacement leading up to the rupture (Bletery and Nocquet, 2023). Our results confirm that, as discussed on informal platforms (Bradley and Hubbard, 2023a,b), the signal is not robust to common-mode filtering. Though this result raises potential concerns on the tectonic origin of the proposed precursory signal, synthetic tests indicate that the common-mode filtering procedure may inadvertently remove an existing signal. Moreover, the collective outcomes of a series of tests we conducted consistently indicate that the signal points to the time, location and slip direction of the impending earthquakes with a statistical significance making very unlikely that the signal solely arises from common-mode noise. The alternative explanation of co-seismic contamination also appears unlikely given that the signal does not appear to be correlated with the co-seismic offsets. Overall, it is difficult to definitely conclude on the origin of the signal. Nevertheless, the interpretation

of the signal as indicative of precursory slip acceleration (Bletery and Nocquet, 2023) remains entirely plausible. Given the potential implications, we encourage others to pursue the investigation in a collaborative effort to confirm or refute the existence of a precursory phase of slow slip leading up to large earthquakes. In that spirit, we are making all our scripts and data available online (see Data and code availability section) for anyone interested to join the effort.

## Acknowledgements

We thank Judith Hubbard, Kyle Bradley, Geoffrey Blewitt, Jeffrey Freymueller, Roland Bürgmann, Jean-Philippe Avouac, Anne Socquet, Yuji Itoh, Kristel Chardard, Martin Vallée, Isabelle Manighetti, Jean-Paul Ampuero, Cédric Twardzik, and François Passelègue for helpful discussion. This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement 949221) and from the French National Research Agency (grant agreement ANR-19-CE310003).

## Data and code availability

All the scripts and data necessary to reproduce the figures presented in this study are available online at the following addresses: https://doi.org/10.5281/zenodo.8064086 (scripts and data from Bletery and Nocquet, 2023), https://doi.org/10.5281/zenodo.11371894 (additional scripts generated for this study), and https://doi.org/10.5281/zenodo.14191538 (scripts corresponding to section 8.2.2, added during the review process).

## Competing interests

The authors declare that they have no competing interests.

## References

* Acosta, M., Passelègue, F. X., Schubnel, A., Madariaga, R., and Violay, M. Can precursory moment release scale with earthquake magnitude? A view from the laboratory. *Geophysical Research Letters*, 46(22):12927-12937, 2019. doi: 10.1029/2019GL084744.

* Amiri-Simkooei, A. Noise in multivariate GPS position time-series. *Journal of Geodesy*, 83:175-187, 2009. doi: 10.1007/s00190-008-0251-8.

* Beaucé, E., Poli, P., Waldhauser, F., Holtzman, B., and Scholz, C. Enhanced tidal sensitivity of seismicity before the 2019 magnitude 7.1 Ridgecrest, California earthquake. *Geophysical Research Letters*, 50(14):e2023GL104375, 2023. doi: 10.1029/2023GL104375.

* Bedford, J. R., Moreno, M., Deng, Z., Oncken, O., Schurr, B., John, T., Báez, J. C., and Bevis, M. Months-long thousand-kilometre-scale wobbling before great subduction earthquakes. *Nature*, 580(7805):628-635, 2020. doi: 10.1038/s41586-020-2212-1.

* Behr, W. M. and Bürgmann, R. What's down there? The structures, materials and environment of deep-seated slow slip and tremor. *Philosophical Transactions of the Royal Society A*, 379 (2193):20200218, 2021. doi: 10.1098/rsta.2020.0218.
Bletery, Q. and Nocquet, J.-M. Slip bursts during coalescence of slow slip events in Cascadia. *Nature communications*, 11(1):1–6, 2020. doi: 10.1038/s41467-020-15494-4.

Bletery, Q. and Nocquet, J.-M. The precursory phase of large earthquakes. *Science*, 381(6655):297–301, 2023. doi: 10.1126/science.adg256.

Blewitt, G., Hammond, W. C., and Kreemer, C. Harnessing the GPS data explosion for interdisciplinary science. *Eos*, 99(10.2019): 485, 2018. doi: 10.1029/2018eo104623.

Bouchon, M., Karabulut, H., Aktar, M., Özalaybey, S., Schmittbuhl, J., and Bouin, M.-P. Extended nucleation of the 1999 M w 7.6 Izmit earthquake. *science*, 331(6019):877–880, 2011. doi: 10.1126/science.1197341.

Bouchon, M., Durand, V., Marsan, D., Karabulut, H., and Schmittbuhl, J. The long precursory phase of most large interplate earthquakes. *Nature geoscience*, 6(4):299–302, 2013. doi: 10.1038/ngeo1770.

Bouchon, M., Marsan, D., Durand, V., Campillo, M., Perfettini, H., Madariaga, R., and Gardonio, B. Potential slab deformation and plunge prior to the Tohoku, Iquique and Maule earthquakes. *Nature Geoscience*, 9(5):380–383, 2016. doi: 10.1038/ngeo2701.

Bradley, K. and Hubbard, J. Earthquake precursors? Not so fast. *Earthquake Insights*, 2023a. doi: 10.62481/310cc439.

Bradley, K. and Hubbard, J. Update on apparent GPS detection of earthquake precursors. *Earthquake Insights*, 2023b. doi: 10.62481/479c2ea4.

Brodsky, E. E. and Lay, T. Recognizing foreshocks from the 1 April 2014 Chile earthquake. *Science*, 344(6185):700–702, 2014. doi: 10.1126/science.1255202.

Bürgmann, R. Reliable earthquake precursors? *Science*, 381 (6655):266–267, 2023. doi: 10.1126/science.adi803.

Caballero, E., Chounet, A., Duputel, Z., Jara, J., Twardzik, C., and Jolivet, R. Seismic and aseismic fault slip during the initiation phase of the 2017 MW= 6.9 Valparaíso earthquake. *Geophysical Research Letters*, 48(6):e2020GL091916, 2021. doi: 10.1029/2020GL091916.

Chanard, K., Fleitout, L., Calais, E., Rebischung, P., and Avouac, J.-P. Toward a global horizontal and vertical elastic load deformation model derived from GRACE and GNSS station position time series. *Journal of Geophysical Research: Solid Earth*, 123 (4):3225–3237, 2018. doi: 10.1002/2017JB015245.

Choi, K., Bilich, A., Larson, K. M., and Axelrad, P. Modified sidereal filtering: Implications for high-rate GPS positioning. *Geophysical Research Letters*, 31(22), 2004. doi: 10.1029/2004GL021621.

Dieterich, J. H. and Kilgore, B. Implications of fault constitutive properties for earthquake prediction. *Proceedings of the National Academy of Sciences*, 93(9):3787–3794, 1996. doi: 10.1073/pnas.93.9.3787.

Dong, D., Fang, P., Bock, Y., Cheng, M., and Miyazaki, S. Anatomy of apparent seasonal variations from GPS-derived site position time series. *Journal of Geophysical Research: Solid Earth*, 107 (B4):ETG-9, 2002. doi: 10.1029/2001JB000573.

Dong, D., Fang, P., Bock, Y., Webb, F., Prawirodirdjo, L., Kedar, S., and Jamason, P. Spatiotemporal filtering using principal component analysis and Karhunen-Loeve expansion approaches for regional GPS network analysis. *Journal of geophysical research: solid earth*, 111(B3), 2006. doi: 10.1029/2005JB003806.

Ellsworth, W. L. and Bulut, F. Nucleation of the 1999 Izmit earthquake by a triggered cascade of foreshocks. *Nature Geoscience*, 11(7):531–535, 2018. doi: 10.1038/s41561-018-0145-1.

Geller, R. J. Earthquake prediction: a critical review. *Geophysical Journal International*, 131(3):425–450, 1997. doi: 10.1111/j.1365-246X.1997.tb06588.x.

Gobron, K., Rebischung, P., Chanard, K., and Altamimi, Z. Anatomy of the spatiotemporally correlated noise in GNSS station position time series. *Journal of Geodesy*, 98(5):34, 2024. doi: 10.1007/s00190-024-01848-z.

Gomberg, J., 2007, C., and Group, B. W. Slow-slip phenomena in Cascadia from 2007 and beyond: A review. *Bulletin*, 122(7-8): 963–978, 2010. doi: 10.1130/B30287.1.

Gualandi, A., Serpelloni, E., and Belardinelli, M. E. Blind source separation problem in GPS time series. *Journal of Geodesy*, 90 (4):323–341, 2016. doi: 10.1007/s00190-015-0875-4.

Helmstetter, A. and Sornette, D. Foreshocks explained by cascades of triggered seismicity. *Journal of Geophysical Research: Solid Earth*, 108(B10), 2003. doi: 10.1029/2003JB002409.

Hirose, H., Kato, A., and Kimura, T. Did short-term preseismic crustal deformation precede the 2011 great Tohoku-oki earthquake? An examination of stacked tilt records. *Geophysical Research Letters*, 51(12):e2024GL109384, 2024. doi: 10.1029/2024GL109384.

Hulbert, C., Rouet-Leduc, B., Johnson, P. A., Ren, C. X., Rivière, J., Bolton, D. C., and Marone, C. Similarity of fast and slow earthquakes illuminated by machine learning. *Nature Geoscience*, 12 (1):69–74, 2019. doi: 10.1038/s41561-018-0272-8.

Jones, L. M. and Molnar, P. Some characteristics of foreshocks and their possible relationship to earthquake prediction and premonitory slip on faults. *Journal of Geophysical Research: Solid Earth*, 84(B7):3596–3608, 1979. doi: 10.1029/JB084iB07p03596.

Kagan, Y. Y. Are earthquakes predictable? *Geophysical Journal International*, 131(3):505–525, 1997. doi: 10.1111/j.1365-246X.1997.tb06595.x.

Kaneko, Y., Nielsen, S. B., and Carpenter, B. M. The onset of laboratory earthquakes explained by nucleating rupture on a rate-and-state fault. *Journal of Geophysical Research: Solid Earth*, 121(8):6071–6091, 2016. doi: 10.1002/2016JB013143.

Kato, A., Obara, K., Igarashi, T., Tsuruoka, H., Nakagawa, S., and Hirata, N. Propagation of slow slip leading up to the 2011 M w 9.0 Tohoku-Oki earthquake. *Science*, 335(6069):705–708, 2012. doi: 10.1126/science.1215141.

Kreemer, C. and Blewitt, G. Robust estimation of spatially varying common-mode components in GPS time-series. *Journal of geodesy*, 95(1):13, 2021. doi: 10.1007/s00190-020-01466-5.

Larson, K. M., Bilich, A., and Axelrad, P. Improving the precision of high-rate GPS. *Journal of Geophysical Research: Solid Earth*, 112 (B5), 2007. doi: 10.1029/2006JB004367.

Latour, S., Schubnel, A., Nielsen, S., Madariaga, R., and Vinciguerra, S. Characterization of nucleation during laboratory earthquakes. *Geophysical Research Letters*, 40(19):5064–5069, 2013. doi: 10.1002/grl.50974.

Lebihain, M., Roch, T., Violay, M., and Molinari, J.-F. Earthquake nucleation along faults with heterogeneous weakening rate. *Geophysical Research Letters*, 48(21):e2021GL094901, 2021. doi: 10.1029/2021GL094901.

Mansinha, L. and Smylie, D. The displacement fields of inclined faults. *Bulletin of the Seismological Society of America*, 61(5): 1433–1440, 1971. doi: 10.1785/BSSA0610051433.

Mao, A., Harrison, C. G., and Dixon, T. H. Noise in GPS coordinate time series. *Journal of Geophysical Research: Solid Earth*, 104 (B2):2797–2816, 1999. doi: 10.1029/1998JB900033.

Martínez-Garzón, P. and Poli, P. Cascade and pre-slip models oversimplify the complexity of earthquake preparation in nature. *Communications Earth & Environment*, 5(1):120, 2024. doi: 10.1038/s43247-024-01285-y.

Mavrommatis, A. P., Segall, P., and Johnson, K. M. A decadal-scale deformation transient prior to the 2011 Mw 9.0 Tohoku-oki
earthquake. *Geophysical Research Letters*, 41(13):4486-4494, 2014. doi: 10.1002/2014GL060139.

Moutote, L., Marsan, D., Lengliné, O., and Duputel, Z. Rare occurrences of non-cascading foreshock activity in southern California. *Geophysical research letters*, 48(7):e2020GL091757, 2021. doi: 10.1029/2020GL091757.

Obara, K. and Kato, A. Connecting slow earthquakes to huge earthquakes. *Science*, 353(6296):253-257, 2016. doi: 10.1126/science.aaf1512.

Ohnaka, M. and Shen, L.-f. Scaling of the shear rupture process from nucleation to dynamic propagation: Implications of geometric irregularity of the rupturing surfaces. *Journal of Geophysical Research: Solid Earth*, 104(B1):817-844, 1999. doi: 10.1029/1998JB00007.

Passelègue, F. X., Latour, S., Schubnel, A., Nielsen, S., Bhat, H. S., and Madariaga, R. Influence of fault strength on precursory processes during laboratory earthquakes. *Fault zone dynamic processes: Evolution of fault properties during seismic rupture*, pages 229-242, 2017. doi: 10.1002/9781119156895.ch12.

Radiguet, M., Perfettini, H., Cotte, N., Gualandi, A., Valette, B., Kostoglodov, V., Lhomme, T., Walpersdorf, A., Cabral Cano, E., and Campillo, M. Triggering of the 2014 Mw7.3 Papanoa earthquake by a slow slip event in Guerrero, Mexico. *Nature Geoscience*, 9 (11):829-833, 2016. doi: 10.1038/ngeo2817.

Rubin, A. M. and Ampuero, J.-P. Earthquake nucleation on (aging) rate and state faults. *Journal of Geophysical Research: Solid Earth*, 110(B11), 2005. doi: 10.1029/2005JB003686.

Ruiz, S., Metois, M., Fuenzalida, A., Ruiz, J., Leyton, F., Grandin, R., Vigny, C., Madariaga, R., and Campos, J. Intense foreshocks and a slow slip event preceded the 2014 Iquique M w 8.1 earthquake. *Science*, 345(6201):1165-1169, 2014. doi: 10.1126/science.1256074.

Ruiz, S., Aden-Antoniow, F., Baez, J., Otarola, C., Potin, B., Del Campo, F., Poli, P., Flores, C., Satriano, C., Leyton, F., et al. Nucleation phase and dynamic inversion of the Mw 6.9 Valparaíso 2017 earthquake in Central Chile. *Geophysical Research Letters*, 44(20):10-290, 2017. doi: 10.1002/2017GL075675.

Scholz, C. H., Sykes, L. R., and Aggarwal, Y. P. Earthquake Prediction: A Physical Basis: Rock dilatancy and water diffusion may explain a large class of phenomena precursory to earthquakes. *Science*, 181(4102):803-810, 1973. doi: 10.1126/science.181.4102.803.

Schurr, B., Asch, G., Hainzl, S., Bedford, J., Hoechner, A., Palo, M., Wang, R., Moreno, M., Bartsch, M., Zhang, Y., et al. Gradual unlocking of plate boundary controlled initiation of the 2014 Iquique earthquake. *Nature*, 512(7514):299-302, 2014. doi: 10.1038/nature13681.

Schwartz, S. Y. and Rokosky, J. M. Slow slip events and seismic tremor at circum-Pacific subduction zones. *Reviews of Geophysics*, 45(3), 2007. doi: 10.1029/2006RG000208.

Socquet, A., Valdes, J. P., Jara, J., Cotton, F., Walpersdorf, A., Cotte, N., Specht, S., Ortega-Culaciati, F., Carrizo, D., and Norabuena, E. An 8 month slow slip event triggers progressive nucleation of the 2014 Chile megathrust. *Geophysical Research Letters*, 44(9): 4046-4053, 2017. doi: 10.1002/2017GL073023.

Tape, C., Holtkamp, S., Silwal, V., Hawthorne, J., Kaneko, Y., Ampuero, J. P., Ji, C., Ruppert, N., Smith, K., and West, M. E. Earthquake nucleation and fault slip complexity in the lower crust of central Alaska. *Nature Geoscience*, 11(7):536-541, 2018. doi: 10.1038/s41561-018-0144-2.

Tian, Y. and Shen, Z.-K. Extracting the regional common-mode component of GPS station position time series from dense continuous network. *Journal of Geophysical Research: Solid Earth*, 121(2):1080-1096, 2016. doi: 10.1002/2015JB012253.

Vallée, M. and Douet, V. A new database of source time functions (STFs) extracted from the SCARDEC method. *Physics of the Earth and Planetary Interiors*, 257:149-157, 2016. doi: 10.1016/j.pepi.2016.05.012.

van den Ende, M. P. and Ampuero, J.-P. On the statistical significance of foreshock sequences in Southern California. *Geophysical Research Letters*, 47(3):e2019GL086224, 2020. doi: 10.1029/2019GL086224.

Voosen, P. Warning signs detected hours ahead of big earthquakes. *Science*, 2023. doi: 10.1126/science adj8753.

Wallace, L. M. Slow slip events in New Zealand. *Annual Review of Earth and Planetary Sciences*, 48:175-203, 2020. doi: 10.1146/annurev-earth-071719-055104.

Wdowinski, S., Bock, Y., Zhang, J., Fang, P., and Genrich, J. Southern California permanent GPS geodetic array: Spatial filtering of daily positions for estimating coseismic and postseismic displacements induced by the 1992 Landers earthquake. *Journal of Geophysical Research: Solid Earth*, 102(B8):18057-18070, 1997. doi: 10.1029/97JB01378.

Wells, D. L. and Coppersmith, K. J. New empirical relationships among magnitude, rupture length, rupture width, rupture area, and surface displacement. *Bulletin of the seismological Society of America*, 84(4):974-1002, 1994. doi: 10.1785/BSSA0840040974.

Williams, S. D., Bock, Y., Fang, P., Jamason, P., Nikolaidis, R. M., Prawirodirdjo, L., Miller, M., and Johnson, D. J. Error analysis of continuous GPS position time series. *Journal of geophysical research*, 2004. doi: 10.1029/2003JB002741.

Zhang, J., Bock, Y., Johnson, H., Fang, P., Williams, S., Genrich, J., Wdowinski, S., and Behr, J. Southern California Permanent GPS Geodetic Array: Error analysis of daily position estimates and site velocities. *Journal of geophysical research: solid earth*, 102 (B8):18035-18055, 1997. doi: 10.1029/97JB01380.

The article *Do large earthquakes start with a precursory phase of slow slip?© 2024 by Quentin Bletery is licensed under CC BY 4.0.