Foreshock-induced slip transients set mainshock nucleation timing

Experimental set-up

A biaxial apparatus, CrackDyn, was used to perform the experiments (Fig. 1a). The apparatus housed two PMMA plates (40 × 10 × 1 cm and 45 × 10 × 1.8 cm). Experiments were conducted on PMMA rather than natural rocks because it provides three main advantages for scaling laboratory observations to natural fault systems. First, owing to its lower elastic stiffness compared with rocks, the state-evolution slip distance, L, and cohesive zone size, X_c, are smaller in PMMA. As a result, a laboratory-scale PMMA fault is dynamically representative of a much larger natural fault, whereas a laboratory rock fault essentially remains of the same scale as its natural counterpart. Second, the reduced elastic properties and fracture energy, G_c, of PMMA ensure that nucleation lengths, ℓ_c and ℓ_∞, remain smaller than the total fault dimension in the experiment, which is rarely the case in rocks, in which nucleation patches can be comparable with or larger than the laboratory fault. Finally, the birefringent properties of PMMA allow for photoelastic visualization, enabling direct observation of nucleation growth and subsequent rupture propagation using high-speed imaging.

The plates were characterized by a static Young’s modulus, E, of 3 GPa and a Poisson’s ratio, ν, of 0.35. A normal load was applied by three vertical pistons by means of steel sample holders. The pistons were supplied by a Enerpac P141 hydraulic pump. Similarly, a shearing load was applied by a single horizontal piston. This piston was supplied by a Top Industrie PMHP-35-1000 hydraulic pump. The force applied by each piston was recorded by a Scaime K13 load cell at 500 Hz. Local strains were recorded by 13 350-ohm strain-gauge rosettes (Micro-Measurements C5K-06-S5198-350-33F) located 3 mm from the sample–sample interface. The strain gauges recorded at 2 MHz and were amplified by a factor of 10 by Elsys SGA-2 MK2 amplifiers. Ten PHILTEC model D100-E2H2PQTS displacement sensors were fixed across the sample–sample interface and recorded the local displacement with a recording frequency of 2 MHz. Thirteen Brüel & Kjær Type 8309 accelerometers were glued horizontally approximately 30 mm from the sample–sample interface and recorded at 2 MHz. These signals were amplified by a NEXUS conditioning amplifier 2692. Finally, an EFFILUX EFFI-LINE3-WTR-600-000-POL-PWR-C light source was used to shine linearly polarized light through the sample. This light then passed through a second linear polarizer before reaching a Phantom TMX 6410 high-speed camera, which recorded 1,280 × 32 pixels at 1 MHz. Considering that PMMA is a birefringent material, this allowed for the use of photoelasticity to track changes in stress across the sample–sample interface^27,37,38,45. A piezoelectric sensor attached to one of the PMMA samples was used as a trigger for an oscilloscope (Picoscope 4224A), which generated a TLL-like signal that triggered the camera and allowed for the synchronization of the other acquisition systems.

Experimental approach

Experiments were begun by applying 100, 150, 200, 250 or 300 bar nominal normal stress (that is, the pressure indicated on the analogue gauge of the pump supplying the vertical pistons). Next, the shear stress was increased by setting a constant flow rate, 3 cm³ min⁻¹, on the pump supplying the horizontal piston. These conditions were kept constant until enough slip had been accumulated such that the displacement sensors were out of range (approximately 0.5 mm). A 1.3 × 0.3 × 10-cm stopper was placed between the lower-most PMMA sample and the horizontal piston. By adjusting this stopper before the experiment such that it was in either a raised or a lowered position, the local loading conditions could be altered such that a wider variety of events could be produced. Each of the five normal stresses were tested with both stopper positions, resulting in ten total experiments and 94 dynamic events.

Data treatment

As the strain gauges were set up in a quarter-bridge configuration, the strain, ε, could be found for the ith strain gauge as,

in which U_i is the voltage reading of the individual strain gauge, U_ex is the excitation voltage, J_amp is the amplification gain and J_f is the gauge factor of the strain gauge.

In this three-strain-gauge rosette, the strain gauges were oriented at 45° from one another, such that the principal strains were then found as

in which ε_v is the vertical direction and ε₃₁₅ and ε₄₅ are oriented 315° and 45° clockwise from vertical (Fig. 1a). Here x refers to the horizontal direction and y refers to the vertical direction.

These strains were then converted to stress, σ, considering Hooke’s law in plane stress conditions,

The values of stress were filtered with a low-pass filter at 10 kHz.

The camera recorded greyscale images in a 1,280 × 32-pixel matrix in which the values ranged from zero (black) to a maximum value (white) of light intensity. Taking the first 20 frames as reference, the future frames were compared with this reference to detect changes in stress during rupture propagation, made possible by the stress-induced birefringence of PMMA and the use of linearly polarized light^37,38. One of the 32 rows (corresponding to the 1,280 pixels closest to the sample–sample interface) was isolated and used for this analysis. The readings of the accelerometers were converted to accelerations using the individually calibrated sensitivities of each sensor. Finally, the readings of the displacement sensors were converted from voltage to gap distance using the calibration provided by PHILTEC. They were then filtered with a low-pass filter at 10 kHz. The nucleation duration was then defined as the interval between the foreshock recorded by the accelerometers and the transition to dynamic rupture, marked by rupture velocities exceeding 10 m s⁻¹ as recorded by high-speed imaging.

Foreshocks and nucleation length scaling

In the present experimental dataset, the onset of nucleation is systematically associated with a foreshock. This feature is specific to this series of experiments and probably reflects the heterogeneous stress distribution along the fault, which varies between PMMA plates. In other laboratory studies, nucleation has been observed to occur without any identifiable foreshock. In such cases, the influence of foreshocks on the nucleation process cannot be investigated and the nucleation size is, to a first order, defined by its theoretical asymptotic value, ℓ_∞ (ref. ²⁰), for understressed faults. The lack of foreshock initiating the quasi-static nucleation phase in other bare-rock-surface experiments can probably be explained by the scaling parameter \(\frac{\Delta T}{{\mu }^{{\prime} }L}\), in which ΔT ≈ μ′δ_a, such that the scaling parameter reads \(\frac{{\delta }_{{\rm{a}}}}{L}\). Because L in PMMA (about 0.1 μm) is one to two orders of magnitude smaller than typical values in bare-surface rock samples (1–10 μm (ref. ⁴⁶)), achieving a comparable influence of a foreshock in rock would require foreshock slips that are 10–100 times larger, or on the order of tens to hundreds of microns. As shown, for example, in ref. ¹⁴, such slip amplitudes are in many cases similar to those of the mainshock, suggesting that the foreshocks observed in those samples are probably much smaller and therefore contribute less to the nucleation dynamics. However, on larger natural faults, the absolute slip accumulated during foreshocks can be substantially greater, potentially enabling foreshocks to exert a dynamically notable perturbation to the EoM. Note, however, that the EoM still describes cases with vanishingly small foreshocks. The limiting case of ΔT → 0 results in ℓ_c → ℓ_∞ for the ageing law (Extended Data Fig. 8). Furthermore, it should be noted that ΔT need not be supplied by a seismic foreshock.

The coexistence of foreshocks and a larger nucleation front can be rationalized by considering the scale dependence of nucleation. At the scale of several micron-sized asperities, local stress concentrations lead to much higher effective normal stresses than the macroscopically applied value. Because the theoretical nucleation size scales inversely with normal stress²⁰,

in which ℓ_b is the elasto-frictional length scale, these local stress amplifications reduce the critical nucleation length at the asperity scale, enabling small slip instabilities that manifest as foreshocks. At the same time, the macroscopic nucleation process spans a much larger portion of the fault, controlled by the bulk stress state and the effective average σ, and can therefore extend over dimensions much greater than the asperity scale. This multiscale interplay explains how foreshocks, driven by localized asperity failure, can coexist with a larger quasi-static nucleation front that evolves towards dynamic rupture, as recently highlighted in the homogenized regime²⁶.

Transient sliding velocity

Although the onset of stable slip is systematically triggered by a foreshock, the associated sliding velocity is too small to be found as an instantaneous measurement. Instead, slip is seen to increase exponentially with time (Extended Data Fig. 2) and can fit with an exponential function. The derivative of this function is found just after the foreshock to provide the transient sliding velocity. A similar procedure is applied to the data digitized for natural earthquakes (Extended Data Fig. 10 and Extended Data Table 1). Note that, in the case of data for natural earthquakes, these data rely on either inversion of the fault slip based on geodetic data or the supposition that the slip during repeater earthquakes is representative of the aseismic slip of the patch⁴⁷. Both methods are subject to several assumptions and potential errors. Illapel 2015 exhibited repeater earthquakes beginning approximately 140 days before the mainshock⁹. The cumulative slip of these repeaters was used to infer aseismic slip, assumed to be representative of quasi-static slip associated with nucleation. Iquique 2014 exhibited repeater earthquakes starting around 270 days before the mainshock⁴. These repeaters were used to infer aseismic slip, which was assumed here to be associated with the nucleation of the mainshock. Alternatively, we could use GNSS data³, which showed substantial movement beginning approximately 17 days before the mainshock to invert for aseismic fault slip. Note that a M_w 6.7 earthquake occurred at the onset of this aseismic slip. However, a time series of this data is not available, such that a calculation of V_min is not possible. Papanoa 2014 exhibited a slow-slip event that began two months before a mainshock⁴⁰. Although GNSS data can be used to infer the slip velocity at the onset of this period, the notion of V_min is unclear, as no known foreshock is associated with the start of the event. Ibaraki-Oki 2008 exhibited repeater earthquakes in the approximately 4 days leading up to the mainshock, with this sequence accelerating substantially approximately 12 h before the mainshock¹⁰. These repeater events were used to estimate aseismic slip assumed to be associated with nucleation. Tohoku-Oki 2011 was preceded by repeating earthquakes, beginning approximately 23 days before the mainshock². Then, a magnitude 7.3 foreshock occurred approximately 2 days before the mainshock. The repeater earthquakes were used to infer aseismic slip. It is considered that the first foreshocks initiated nucleation, resulting in a slip velocity peak followed by a slip velocity minimum (Extended Data Fig. 11). Next, the M7.3 probably had a strong influence on the nucleation dynamics, resulting in another slip velocity peak and a new slip velocity minimum (Extended Data Fig. 11). It is possible that the initial foreshocks are linked to the nucleation of the M7.3, which itself governs the nucleation of the M9.0; however, note that other authors consider that the associated pre-slip occurred away from the mainshock rupture zone⁴⁸. Here the initial foreshocks are considered to initiate the nucleation of the M7.3 event, which itself is considered to initiate the nucleation of the M9.0 mainshock. İzmit 1999 exhibited a 44-min-long foreshock sequence^1,11. However, the exact nature of this sequence is controversial and it is unclear whether there was a precursory aseismic driving process. Here the assumption is made that these are repeating earthquakes representative of aseismic slip, but it should be emphasized that this may not be the case. Valparaíso 2017 was observed to be preceded by repeater-type seismicity and trenchward movement observed in GNSS data beginning 2 days before the mainshock⁶. Although fault slip can be inverted from the GNSS data, the full time series is not available. Instead, a simple average of fault slip over the time interval was used to provide an upper-bound estimate of V_min. Yangbi 2021 was associated with migrating foreshocks and repeating earthquakes, with a rapid migration beginning 30 min before the mainshock¹³. During this period, GNSS data allows for the inversion of aseismic slip, assumed to be associated with nucleation. Unfortunately, the time series of this data is not available, so a simple average is again used to provide an upper-bound estimate of V_min.

Foreshock magnitude

The foreshocks at the origin of all observable nucleation phases in these experiments are located using error minimization based on arrival times. This localization allows for the assessment of the event magnitude by means of an averaged acceleration, A_fs, found as⁴⁹,

in which N_a is the number of accelerometers, r is the distance from the sensor to the event location (normalized by a reference sphere of 0.01 m) and \({A}_{i}^{\text{max}}\) is the maximum acceleration measured by the sensor during the event.

Unified framework for earthquake nucleation from laboratory to nature

The EoM of the crack originates from the Griffith’s energy balance, which together with a set of simplifying assumptions, can allow for the solving of the crack front propagation⁵⁰. Recently, the EoM framework has been applied to crack-like slip transients on rate-and-state frictional faults^20,22,51. Here we follow the EoM treatment by Garagash²², who considers a plane-strain or anti-plane shear crack of half-length ℓ, propagating at instantaneous velocity ν_r along a planar 1D fault, characterized by the ambient sliding velocity V₀ and ‘overstress’, the difference between the initial stress (ratio) f₀ = τ₀/σ and the steady-state level of friction at the ambient sliding velocity, Δf₀ = f₀ − f_ss(V₀). The fault is embedded in an infinite elastic medium and its rate-and-state-dependent resistance to sliding³² is governed by three parameters: the direct effect coefficient a, the state-evolution coefficient b and the characteristic slip distance L over which the state of the sliding surface is renewed (not to be confused with the apparent slip-weakening distance of the rupture, δ_c). Evolution of the state is often modelled in the framework of either ageing¹⁶ or slip⁵² laws, which manifest in different apparent weakening behaviours of strength with incurred slip in a transient (for example, behind a propagating rupture front). Specifically, a transient governed by the ageing law corresponds to a rate of weakening with slip (\(W=-\frac{{\rm{d}}\tau }{{\rm{d}}\delta }\)), which is approximately independent of the strength of the transient^22,53, W ≈ bσ/L, whereas for the slip law, the weakening rate increases with the transient strength, W ≈ Δf_pσ/L, in which Δf_p is the maximum departure of friction from the steady state during the transient (that is, ‘friction breakdown’). Conversely, the apparent slip-weakening distance δ_c required to ‘break’ the strength by amount Δf_pσ from the peak to the steady-state ‘residual’ value at the above values of the weakening rate W is an increasing function of Δf_p, \({\delta }_{{\rm{c}}}=\left(\frac{\Delta {f}_{{\rm{p}}}}{2b}\right)L\), for the ageing law, whereas it is invariant, δ_c ≈ L, for the slip law. Considering that friction experiments have shown the relevance of the ageing law for PMMA⁵⁴, we use the ageing law for the rest of the analysis. However, previous studies have indicated that the slip law may be more appropriate for natural fault material⁵⁵, conversely to PMMA. Note that the model presents a qualitatively similar behaviour using the slip law (Extended Data Fig. 6). For consistency and simplicity, we extend our results to natural earthquakes using the ageing law rather than the slip law (Fig. 4d).

EoM framework

The EoM considers the balance between the energy release rate \(G={\overline{K}}^{2}/(2\bar{\mu })\) into the propagating crack tip with the fracture energy of the frictional breakdown process G_c. Here \(\overline{K}=k({\nu }_{{\rm{r}}})K\) is the dynamic SIF obtained from the static value K multiplied by the wave-mediated dynamic prefactor k(ν_r) and \(\overline{\mu }=g({\nu }_{{\rm{r}}}){\mu }^{{\prime} }\) is the apparent dynamic shear modulus formed from the static value μ′ = μ and dynamic prefactor g(ν_r). These dynamic factors are mostly inconsequential during the nucleation phase characterized by ν_r ≪ c_s, with k ≈ g ≈ 1, and full expressions²² are omitted here for simplicity. The fracture energy of a rate-and-state fault can be defined as the breakdown work G_c = ∫Δfσdδ expended near the rupture front in dissipating the strength excess over the steady-state value, Δf = f − f_ss(V), from its peak value Δf_p, attained at the front, to zero some distance behind the front over some slip distance δ relatable to the state-evolution slip distance L. Because the peak friction, which is approximately alnν_r, and the peak friction breakdown, Δf_p ≈ blnν_r, are increasing with increasing rupture velocity by the direct and state-evolution effects, respectively, so is the fracture energy G_c. Contrary to the energy release rate G, which materially depends on ν_r only in the seismic range, the dependence of the fracture energy G_c on ν_r is persistent over the entire range of rupture speeds, from quasi-static to dynamic. For the ageing law used here, Garagash²² gives

in which Li₂ is the dilog function. This expression further simplifies to \({G}_{{\rm{c}}}\approx (1/2){\left(\frac{\Delta {f}_{{\rm{p}}}}{b}\right)}^{2}b\sigma L\) for transients with Δf_p/b ≫ 1 . The dependence on rupture speed originates from that of the peak friction breakdown Δf_p(ν_r), given implicitly by equation (2.17) in ref. ²². Here we use an explicit expression,

in which \({\bar{\nu }}_{0}=\exp \left(-\frac{\Delta {f}_{0}}{b}\right)\frac{{\mu }^{{\prime} }}{b\sigma }{V}_{0}\) is the characteristic rupture velocity embodying initial fault conditions (ambient sliding velocity V₀ and initial overstress Δf₀ = f₀ − f_ss(V₀)), κ₀ ≈ 0.838 and ProductLog₋₁(x) denotes the −1 branch of the Lambert productlog function providing the real-valued solution of neⁿ = m for −1/e < m < 0. The above explicit inversion of Garagash’s equation (2.17) for Δf_p is exact for the slip law case and is an excellent approximation for the ageing law case considered here. Indeed, Garagash²² observes that Δf_p(ν_r) is essentially indistinguishable for the two laws when plotted in his Fig. 4a.

Reasonable (zeroth-order) approximations, \(\Delta {f}_{{\rm{p}}}\approx b\mathrm{ln}({\nu }_{{\rm{r}}}/{\bar{\nu }}_{0})\) and \({G}_{{\rm{c}}}=(1/2){(\mathrm{ln}({\nu }_{{\rm{r}}}/{\bar{\nu }}_{0}))}^{2}b\sigma L\), allow us to readily glance at a simplified (logarithm-based) relation of the peak friction breakdown and the associated fracture energy on the rupture speed across the entire range of slip transients, from slow to fast.

The crack energy balance G = G_c can be further recast as K = K_c, in which the static SIF K matches the effective fault toughness \({K}_{{\rm{c}}}=(\sqrt{g}/k)\sqrt{2{\mu }^{{\prime} }{G}_{{\rm{c}}}}\), with the prefactor \((\sqrt{g}/k)({\nu }_{{\rm{r}}})\) encapsulating the wave-mediated effects and G_c(ν_r) the quasi-static-to-dynamic range of the fracture energy.

The SIF K = K_bg + K_fs consists of the contributions from the background stress, K_bg, and the foreshock, K_fs, which can be approximated by the solutions for a bi-wing crack of half-length ℓ loaded by an effective uniformly distributed stress drop Δτ_eff and a localized, hypocentral point force ΔT, respectively,

We note that an edge-crack geometry and corresponding slightly modified expressions for the SIF components can be alternatively used to describe the experimental ruptures, which start proximally to one of the experimental fault ends. However, such a modification is not expected to greatly improve the prediction of an already much simplified model. To close-form the EoM, we now inform the SIFs with expressions for the effective stress drop Δτ_eff and hypocentral force ΔT. The hypocentral (foreshock) force ΔT, by a dimensional argument, is

in which δ_a is the foreshock ‘asperity’ slip and C ≈ 0.3 is a numerical coefficient, which will be used as a fitting parameter. The effective (uniform) stress drop Δτ_eff is the energetically equivalent measure of the actual, spatially varying stress drop along the crack defined as,

the drop from the initial stress f₀σ to the ‘residual’ stress f_ss(V)σ set by the frictional steady-state

as parametrized here in terms of the fault initial state: f₀, V₀ and overstress Δf₀ = f₀ − f_ss(V₀). The spatial distribution of sliding velocity V = V(x) along the crack then defines the spatial distribution of the stress drop. The energetically equivalent Δτ_eff can be found from matching its SIF contribution (that is, K_bg) to that of Δτ(V(x)), that is,

Alternatively, we rewrite

in which the effective sliding velocity V_eff is implicitly defined by the second equality when the sliding velocity V(x) is known. Approximating the latter from the solution for the slip on the propagating crack loaded by Δτ_eff and ΔT (Section 3(b) of ref. ²²) and solving for V_eff ultimately results in

Finally, the EoM

together with the above expressions for K_c(ν_r) and V_eff(ν_r, ℓ) can be solved for ν_r as a function of crack length ℓ. The relation between these two parameters can be integrated to solve for time. Evolution of the fracture toughness, fracture energy and sliding velocity also result.

In the main text, we simplify equation (16) by introducing \({K}_{{\rm{bg}}}=\Delta \tau ({V}_{{\rm{eff}}})\sqrt{\pi {\ell }}\). K_bg is increasing with crack half-length ℓ and scales with the rate-dependent transient stress drop Δτ—a function of the initial fault overstress parameter Δf₀ defined as the distance between the background stress ratio f₀ from the steady-state friction value f_ss(V₀) at the ambient sliding velocity V₀—and transient rupture velocity ν_r. In our experiments, the fault is probably neutrally or mildly understressed (Δf₀ ≱ 0; Fig. 3c and Extended Data Fig. 12). Furthermore, we introduce \({K}_{{\rm{fs}}}=\Delta T/\sqrt{\pi {\ell }}\), which denotes the SIF contribution from the foreshock. It decreases with the crack half-length ℓ and scales with a localized hypocentral Coulomb force of magnitude ΔT ∝ μ′δ_a, generated by the impulse associated with the foreshock slip, δ_a, that initiates the nucleation process.

We note that the classical nucleation crack half-length ℓ_∞ (ref. ²⁰) (equation (5)) emerges from the EoM in a formal asymptotic limit of diverging, but not yet seismic, sliding velocity and rupture speed. Indeed, in this case, the SIF, dominated by the stress drop contribution, \(K\approx {K}_{{\rm{bg}}}\approx (b-a)\sigma \sqrt{\pi {\ell }}\times \mathrm{ln}({V}_{{\rm{eff}}}/{V}_{0})\), balances the asymptotic expression for the fault toughness \({K}_{{\rm{c}}}\approx \sqrt{{\mu }^{{\prime} }b\sigma L}\times \mathrm{ln}({\nu }_{{\rm{r}}}/{\bar{\nu }}_{0})\). Given the proportionality of the sliding velocity and rupture speed (equation (15)), the two logarithmically diverging terms in the above are identical to the leading order, resulting in \((b-a)\sigma \sqrt{\pi {\ell }}\approx \sqrt{{\mu }^{{\prime} }b\sigma L}\), with the solution ℓ = ℓ_∞. Conditions under which the theoretical nucleation length ℓ_∞ describes the actual nucleation process depend on the validity of the assumption of an unbounded sliding velocity at nucleation, which allows for the neglecting of finite contributions of the foreshock ΔT and understress/overstress Δf₀ to the SIF, as well as the neglecting of the difference in the scaling factors under the logarithms in expressions for K_bg and K_c in the above asymptotic analysis. In reality, these various contributions would have to be evaluated and compared at a finite rupture speed and sliding velocity of the nucleation process (a fraction of seismic values) to validate the applicability of the asymptotic solution.

Illustrative examples of EoM solution for the sliding velocity histories resembling the experimentally observed, foreshock-mediated rupture nucleation regimes are shown in Extended Data Fig. 6 and Fig. 3b using the parametrization explained in the following. See Extended Data Table 2 for a tabulation of the main variables.

EoM solution in parametric space

We use a suite of several simulations of the EoM to investigate the dependence of the foreshock-mediated nucleation process: foreshock, transient slow slip, acceleration to mainshock, as encapsulated by the nucleation time Δt_c, nucleation length ℓ_c and the transient minimum sliding velocity V_eff,min, on the foreshock force ΔT, initial fault conditions (Δf₀ and V₀) and a set of fault elasto-frictional constitutive parameters. We use a normalized EoM formulation that allows for the reduction of the parametric space to a non-dimensional foreshock force, \(\frac{\Delta T}{{\mu }^{{\prime} }L}\), initial overstress, \(\frac{\Delta {f}_{0}}{b}\), ambient sliding velocity, \(\frac{{V}_{0}}{(b\sigma /{\mu }^{{\prime} }){c}_{{\rm{s}}}}\), and a single constitute parameter, the ratio of the direct to state coefficients, a/b. We present results for the PMMA-like value a/b ≈ 0.55 (ref. ⁵⁶) while also performing a limited suite of solutions for the rock-like value a/b ≈  0.75 (ref. ³⁵) to further validate the emerging correlations discussed below.

Nucleation length, nucleation time, transient minimum velocity and inversion of fault-state-evolution distance

The EoM suite of solutions for the normalized nucleation length versus nucleation time is shown in Extended Data Fig. 8 for a wide variety of the non-dimensional loading parameters foreshock force \(\frac{\Delta T}{{\mu }^{{\prime} }L}\) (implicit), initial overstress \(\frac{\Delta {f}_{0}}{b}\) and nondimensional ambient sliding velocity \(\frac{{V}_{0}}{(b\sigma /{\mu }^{{\prime} }){c}_{{\rm{s}}}}\). Length and time to nucleation decrease with both increasing foreshock impulse and increasing fault initial overstress. The theoretical nucleation length ℓ_∞ (refs. ^20,57) serves as the upper bound, which is approached asymptotically on understressed (Δf₀ < 0) faults when the foreshock force decreases to the Δf₀-dependent threshold value for the transient arrest (and as time to nucleation diverges).

Performed EoM simulations show a collapse of time-to-nucleation Δt_c versus minimum effective sliding velocity, V_eff,min, for a wide variety of the non-dimensional loading parameters, \(\frac{\Delta T}{{\mu }^{{\prime} }L}\) and \(\frac{\Delta {f}_{0}}{b}\), and ambient sliding velocity \(\frac{{V}_{0}}{(b\sigma /{\mu }^{{\prime} }){c}_{{\rm{s}}}}\) (Extended Data Fig. 9a). This collapse manifests two asymptotic behaviours in which Δt_c ≈ L/V₀ is a constant for transients with very small minimum sliding velocity, comparable with the initial fault value, that is, when V_eff,min ≈ V₀, and Δt_c ≈ (L/V_eff,min)ln(V_eff,min/V₀) when V_eff,min ≫ V₀. Essentially, to the first order, Δt_c is inverse in V_eff,min with the low-velocity cut-off at about V₀. Small departure from the collapsed curve on Extended Data Fig. 9a corresponds to the solutions with appreciable understress Δf₀ < 0. This second-order behaviour is more apparent when we remove the dominant (inverse in V_eff,min) behaviour when plotting the product Δt_cV_eff,min in Extended Data Fig. 9b.

The overall relation between Δt_c and V_eff,min (Extended Data Fig. 9a) can be well modelled by

with the transition function

between the small-velocity, \({\mathcal{D}}\approx {C}_{0}{\mathcal{V}}\), and large-velocity, \({\mathcal{D}}\approx {C}_{1}\mathrm{ln}\left(\frac{{\mathcal{V}}}{{{\mathcal{V}}}_{1}}\right)\), asymptotes, and fitting constants C₀, C₁, \({{\mathcal{V}}}_{1}\) and α found as 0.64174, 1.33245, 0.01332 and 2.15680, respectively. The asymptotes follow by considering the time required to achieve the maximum stable crack length, ℓ_∞, at the minimum rupture velocity, ν_r,min, itself dependent on V_eff,min by equation (15).

To perform the fit of the EoM empirical relation of the time-to-nucleation Δt_c and the average minimum sliding velocity V_eff,min (equation (17)) to the correlation of experimental events Δt_c with the minimum value of the sliding velocity V_min measured at a fixed point on the fault, we need to relate the pointwise value V_min to the crack-average one V_eff,min. To accomplish this, we use the EoM approximate solution for the sliding velocity distribution²² V(x) = (V_eff/2)(1 − (x/ℓ)²)^−1/2 and model the sliding velocity measured at a fixed location on the fault proximal to the foreshock as the hypocentral value of the EoM distribution. In other words, we use the V_eff,min/2 of the model as the proxy for the experimental V_min.

Taking a = 0.016 and b = 0.029 for the PMMA laboratory fault⁵⁶, equation (17) was fit to the laboratory data, yielding a fault ambient sliding velocity V₀ = 1.35 × 10⁻⁶ m s⁻¹ and the state-evolution slip distance of the interface L = 0.19 μm comparable with approximately 0.4 μm inferred from velocity-stepping experiments on a PMMA fault⁵⁶.

Similarly, in Fig. 4d, equation (17) was used to prepare possible fits for the natural earthquake data assuming granite friction³⁵ with a = 0.014 and b = 0.019 and natural fault ambient sliding velocity V₀ = 10⁻¹² m s⁻¹. We note that the transient (minimum) sliding velocity values inferred for the natural earthquakes in our dataset (Fig. 4d) are on the order of or exceed the typical plate convergence rate, that is, V_min ≳ 10⁻⁹ m s⁻¹. Thus, any suitable choice of the ambient sliding velocity assuming an initially ‘locked’ fault, that is, V₀ ≪ 10⁻⁹m s⁻¹, will not have a substantial effect on the EoM fit for the natural events.

Transient minimum velocity versus foreshock force versus foreshock slip

The EoM empirically predicts a collapse of V_eff,min and ΔT only in the case of large overstress (Extended Data Fig. 12). In those cases, we find \(\frac{{V}_{\mathrm{eff},\text{min}}}{{V}_{0}}\approx 10\times {{\rm{e}}}^{\frac{\Delta T}{0.51{\mu }^{{\prime} }L}}\), in which ΔT = Cμ′δ_a. We use this relation together with the fitting constant C = 0.3 to produce the upper-bound dashed curves in Fig. 3c for different values of L bracketing the best-fit value of 0.19 μm. Indeed, the EoM does not indicate that there is a one-to-one relationship between V_eff,min and the hypocentral force, ΔT. Instead, V_eff,min also depends on the overstress. Rather, V_eff,min and its experimental proxy V_min encapsulate the loading conditions on the fault from both prestress and the foreshock (asperity slip) and cannot be reduced to one or the other.

EoM modelling of experimental sliding velocity history

The model parameters used to produce the experiment-reflecting results for transient sliding velocity evolution of Fig. 3b are the same as used above for PMMA with the addition of the shear wave speed taken as 1,345 m s⁻¹, the apparent shear modulus taken as 1.24 GPa and the normal stress taken as 0.8333 MPa (the average of the normal stresses given in Fig. 2). The asperity slip δ_a is converted into an equivalent hypocentral Coulomb force ΔT using a single asperity equivalent point force, ΔT = Cμ′δ_a with C = 0.3.