Foreshock-induced slip transients set mainshock nucleation timing

2026-05-21Fryer,Barnaby,Garagash,Dmitry,Lebihain,Mathias,Passelègue,François

Nature 2026年8115期

关键词：公式

Abstract

Foreshocks are sometimes observed before earthquakes 1 , 2 , 3 , 4 , 5 , 6 , 7 , 8 , 9 , 10 , 11 , 12 , 13 , yet their role in controlling rupture nucleation remains unclear 1 , 11 , 14 . Classical models often assume that nucleation arises from slow, quasi-static slip governed primarily by fault weakening 15 , 16 , 17 , 18 , 19 , 20 , 21 , typically neglecting impulsive precursory events. Here we show, using laboratory experiments and a rate-and-state-based Griffith-like rupture framework 22 , that foreshocks, when they occur at the onset of or during nucleation, can fundamentally regulate earthquake initiation. We find that the slip burst induced by foreshocks imparts a transient sliding velocity, V min , whose magnitude is set by foreshock size and which robustly predicts both nucleation duration and spatial length. Larger foreshocks generate higher V min and trigger a more rapid transition to dynamic rupture, whereas smaller foreshocks produce long-duration quasi-static growth and very small impulses lead to ruptures entirely arresting. Extending our theoretical framework to tectonic faults, we show that foreshock and associated slow-slip sequences preceding natural earthquakes seem to follow the same scaling. These observations allow us to constrain realistic characteristic nucleation slip distances of 0.3–3.0 mm, orders of magnitude smaller than those inferred for dynamic rupture 23 . Our results demonstrate that foreshock-induced transients set the timing and potential detectability of earthquake nucleation 24 .

Main

Understanding the nucleation of frictional ruptures is a central challenge in earthquake mechanics 18 . Beyond earthquake physics, nucleation processes are generic to frictional instabilities. Nucleation refers to the early processes that initiate fault rupture, marking the transition from a stable, locked fault to an unstable, dynamically slipping one. Coseismic rupture is theoretically expected to be preceded by quasi-static slip, growing in space and time until it reaches a critical size that may trigger dynamic failure 15 , 16 , 17 , 18 , 19 , 20 , 21 . This behaviour is observed in both rate-and-state 16 , 20 , 22 , 25 and slip-weakening models 15 , 17 , 19 , 26 . These foundational studies suggest that the nucleation length, a critical measure of the size of the slipping zone before dynamic rupture, is primarily governed by the fracture energy of the fault interface, G c , and the weakening rate within the nucleation region 19 . However, the physical parameters that control the duration and detectability of this stable phase remain poorly understood 25 , 27 . As a result, it is still unclear under which conditions geophysical observations might reveal precursory signals. Some earthquakes seem to be preceded by a prolonged nucleation stage 2 , 3 , 4 , 6 , 8 , 24 , potentially accompanied by foreshocks 1 , 3 , 4 , 6 , 8 , 28 , whereas many others exhibit no detectable dynamic precursor activity at all. These contrasting observations highlight the need for further investigation into the various physical mechanisms that may govern earthquake nucleation. Key open questions include: (1) whether all earthquakes exhibit quasi-static nucleation; (2) how foreshocks affect nucleation; and (3) whether the slip distance governing dynamic rupture also applies to nucleation.

We report laboratory experiments that, for the first time to our knowledge, directly capture how foreshocks control the nucleation of earthquakes. A small initial slip event sets the transient sliding velocity, which in turn dictates the duration and scale of the nucleation phase before dynamic rupture. Our results are supported by a theoretical framework based on a rate-dependent Griffith-like equation of motion (EoM), formulated using rate-and-state friction. Finally, we extend this theoretical framework to natural earthquakes preceded by foreshocks, enabling constraint of the characteristic nucleation parameters of large earthquakes, challenging standard estimates.

Two decades of nucleation duration

We conducted experiments using a biaxial loading apparatus (Fig. 1a; see Methods for details). Three representative nucleation phases from a single experiment are shown in Fig. 1b–e. Despite having nearly identical initial stress conditions (Fig. 1b), the duration of the nucleation phase varied by more than an order of magnitude, ranging from 0.6 to 12.7 ms (Fig. 1c–e), and reached maximum values of about 80 ms across the experimental suite, highlighting the intrinsic complexity of nucleation dynamics.

**Fig. 1: Large variations in nucleation duration despite similar stress states.**

Note that the three selected representative nucleation phases are observed to initiate with a foreshock (Fig. 1c–e and Extended Data Figs. 1g–i and 2a–c). In many cases, further foreshocks can be clearly identified before the onset of the main rupture (Extended Data Fig. 1g,h). The initiation of nucleation with a distinct foreshock is a consistent feature across the experimental dataset. The primary exceptions are events in which no discernible nucleation phase is present; in these instances, the initial foreshock evolves directly into dynamic rupture, without the intermediate stage of detectable quasi-static slip-front propagation (Extended Data Fig. 3).

Weakening rate and nucleation duration

As slip, δ, initiates within the nucleating patch, the shear stress, τ, begins to decrease, reflecting fault weakening. Initially, the slip-weakening rate, fracdtau ddelta, is steep, meaning that a substantial stress drop occurs over a small amount of slip. However, the weakening rate progressively decreases as slip continues, eventually approaching a quasi-steady residual shear stress level (Fig. 2). Previous studies have shown that sharper weakening rates result in smaller nucleation lengths 17 , 19 . These studies generally indicate that weakening rate has a similar effect on nucleation duration, which further depends on, for example, loading rate 29 . Yet, the example events show similar weakening rates despite an order of magnitude difference in nucleation duration (Fig. 2).

**Fig. 2: Comparable weakening rates for an order of magnitude difference in nucleation duration.**

Moreover, the two orders of magnitude variation in nucleation durations observed across the full experimental suite cannot be easily attributed to differences in weakening rate. In fact, the initial weakening seems to be insensitive to pressure, that is, non-frictional, as the drop in strength is essentially constant despite substantial differences in normal stress for many events.

Foreshocks and the onset of sliding

Although the stress criticality does not seem to vary greatly between cases of slow and rapid nucleation (Fig. 1b), the initial impulse imparted to the nucleating patch does (Extended Data Fig. 4). The initial foreshock generates a burst of sliding velocity along the fault. This rapidly decays to a minimum transient sliding velocity, V min (Fig. 3a,b). The value of V min is found to be positively correlated with the size of the foreshock, such that larger foreshocks produce higher transient sliding velocities (Fig. 3b,c). This relationship suggests that the initial energy imparted by the foreshock plays a key role in controlling the evolution and duration of the nucleation phase.

**Fig. 3: Foreshock size predicts the transient sliding velocity of the nucleating patch.**

Effectively, a clear correlation emerges between nucleation duration and foreshock size (Fig. 4a): events with longer nucleation phases are typically preceded by smaller foreshocks and accumulate greater amounts of quasi-static slip before transitioning to dynamic rupture. This trend is not only evident in the present experimental dataset but also seems consistent with observations from natural earthquakes (Fig. 4a). A general pattern of nucleation behaviour can thus be identified (Fig. 3a,b). The sequence begins with a foreshock that triggers a burst in sliding velocity locally along the fault and crack expansion. This sliding velocity then decays to a minimum, that is, the so-called transient sliding velocity V min , corresponding to a quasi-static slip phase, before both sliding velocity and rupture velocity accelerate exponentially (Extended Data Fig. 2), ultimately leading to dynamic rupture onset.

**Fig. 4: Transient sliding velocity predicts nucleation duration and length.**

Although natural earthquake rupture may not always be triggered by discrete foreshocks, larger foreshocks lead to higher transient sliding velocities, V min , both during and immediately following such precursory activity. This trend is clearly evident in our experimental data, in which larger foreshocks correlate with higher values of V min (Fig. 3c).

Equation of motion

We now seek to reconcile our experimental observations with a state-of-the-art rupture model, grounded in rate-and-state friction. Garagash 22 developed an EoM linking rate-and-state parameters and loading conditions to describe the nucleation and propagation of frictional ruptures. It reads:

(1)

Here all terms are implicit functions of the frictional parameters (a, b, L), in which a, b and L and are the direct effect, evolution effect and state-evolution slip distance, respectively. K bg is the stress intensity factor (SIF) owing to the background stress. It depends on the initial overstress Δf 0 at ambient sliding velocity V 0 , the crack half-length ℓ and rupture velocity ν r , following rate-and-state friction. K fs denotes the SIF contribution from the foreshock and depends on ℓ and a localized Coulomb force of magnitude ΔT∝μ′δ a , in which μ′ is the apparent shear modulus and δ a is the asperity slip. Ultimately, the two SIF contributions must balance the fault toughness K c , which increases with ν r , owing to the rate-dependent frictional response. The exact form of the SIFs is not required to qualitatively match experimental observations, as long as K bg increases with ℓ but K fs decreases with ℓ. This defines a ‘well’ in the energy available to drive the rupture at some intermediate crack length 30 , 31 , which, when balanced with K c (ν r ), sets a transient minimum in the rupture propagation and sliding velocities (Methods).

Nucleation dynamics are therefore controlled by two parameters, Δf 0 and ΔT, which scale the aforementioned dependence of K bg and K fs on the crack length, respectively. As shown in Extended Data Fig. 6, the EoM predicts three possible nucleation scenarios: (1) a transient deceleration followed by arrest without the nucleation of dynamic rupture; (2) deceleration to V min followed by a prolonged quasi-static phase during which both slip and rupture velocity accelerate, eventually transitioning to dynamic rupture; or (3) near-immediate acceleration into dynamic rupture, without a quasi-static phase. All three behaviours are observed in our experimental dataset (Extended Data Fig. 7, Fig. 3a and Extended Data Fig. 3) and can be rationalized using the EoM. Indeed, following equation (1), rupture initiation is initially driven by the impulsive force ΔT. As the crack grows, the foreshock contribution K fs decreases with increasing ℓ, promoting rupture deceleration. However, the background contribution K bg increases with ℓ and tends to accelerate the crack. Crack growth, in turn, raises the sliding velocity V, which increases the transient stress drop, owing to the rate-weakening nature of the interface, and feeds K bg . This feedback can sustain or arrest rupture depending on foreshock magnitude, represented by ΔT, as it defines early-stage sliding velocity and thus the transient strength of the interface. When ΔT is too small, the increase in K bg with increasing crack half-length ℓ is insufficient to offset the associated decrease in K fs and rupture arrests (regime 1). By contrast, a large ΔT generates high transient sliding velocities, leading to rapid growth of K bg and immediate transition to dynamic rupture (regime 3). Between these extreme cases, a moderate ΔT induces a slower but self-sustained increase in K bg , reinforced by the positive feedback of ℓ and V, resulting in delayed dynamic rupture (regime 2). Here the overstress Δf 0 primarily sets the baseline value of the transient stress drop and thus controls the threshold values of ΔT that determine the transition between nucleation regimes.

Moreover, the EoM offers an explanation for the observed variability of the nucleation length ℓ c and the lack of a clear correlation with normal stress σ (Fig. 4b), which is otherwise expected in nucleation models 19 , 20 , 32 and reported in some experimental studies 27 , 33 , 34 , neither of which generally consider impulses. The predicted nucleation size of the EoM for ageing frictional faults is bounded by the maximum value ell _infty ~ (1)/(pi )frac1(1-a/b) 2 fracmu ^prime Lbsigma (ref. 20 ). This terminal value is, however, approached by nucleation transients driven by a vanishingly small impulse on a neutrally or slightly understressed fault (Extended Data Fig. 8). Otherwise, both the overstress and the impulse imparted by the foreshock act to materially reduce the nucleation length compared with ℓ ∞ by as much as an order of magnitude (Extended Data Fig. 8).

In summary, the EoM offers a physical explanation for the observed rupture behaviour during nucleation. The advancing crack tip draws energy from both the inherent stress state along the fault and the stress concentration (or redistribution) generated by the foreshock. However, these two contributions are not necessarily readily constrainable in both laboratory and natural settings. Instead, as these two components are the only contributors to the energy balance of the EoM, their combination results in one unique nucleation behaviour or crack development path (Fig. 3b). This unique crack development path is characterized by the minimum transient sliding velocity V min , which then, as a single parameter, dictates nucleation behaviour (Fig. 4c and Extended Data Fig. 9a).

Transient slip rate controls nucleation

The nucleation length of the experimental events is positively correlated with the nucleation duration (Fig. 4b). This trend holds not only in our experiments but also in larger-scale (about 1 m) rock friction experiments without a distinguishable foreshock and for quasi-statically slipping patches of natural earthquakes, in which the onset of nucleation is identified by a foreshock. These foreshocks provide an empirical basis for estimating nucleation duration, Δt c . Larger impulsive forces (that is, larger foreshocks) result in higher transient minimum sliding velocities V min (Fig. 3b,c) and shorter nucleation duration Δt c and nucleation length ℓ c (Fig. 4c). At ‘large’ sliding velocities V min , Δt c scales as the inverse of V min , Delta t_cpropto fracLV min. However, a progressive shift from this inverse power-law behaviour to a constant nucleation duration Δt c ≅ Δt 0 is predicted at the transition between regimes 1 and 2, when V min approaches the ambient sliding velocity V 0 (Extended Data Fig. 9c). The relationship between Δt c and V min is quantitatively predicted by the EoM using the frictional properties of polymethylmethacrylate (PMMA) and realistic laboratory values for the state-evolution slip distance L (equation (17)).

Notably, the observed dependence of nucleation duration on transient sliding velocity seems to extend to quasi-static slip patches preceding certain natural earthquakes (Fig. 4d). The slope of this relationship is consistent between our laboratory experiments and natural events, but the natural earthquake data are offset by several orders of magnitude. This offset implies that, for a given transient velocity, natural faults exhibit much longer nucleation durations than those observed in the laboratory. This shift is attributable to differences in frictional properties, specifically much larger state-evolution slip distances and smaller values of the rate-weakening coefficient b − a for natural faults, which may depend on fault lithology 35 . Furthermore, part of the discrepancy may arise from the uncertainty in estimating transient sliding velocities for natural earthquakes; our values, based on indirect constraints and the assumption of exponentially increasing slip, probably overestimate the V min during nucleation. Indeed, these data contain notable uncertainty because transient sliding velocities for natural earthquakes must be inferred indirectly from either geodetic slip inversions or by assuming small repeater earthquakes recover accrued aseismic slip deficit (Extended Data Fig. 10). Nevertheless, variations in the a/b ratio arising from possible differences in rheology and/or variations of state-evolution slip distance L, possibly related to fault maturity 36 , are sufficient to explain the discrepancies observed among the different earthquakes analysed (Fig. 4d).

Our experimental results provide new insight into the physics of earthquake nucleation and allow us to revisit several key questions about the onset of seismic rupture. (1) Are all earthquakes preceded by a quasi-static nucleation phase? Although most of the dynamic events in these experiments exhibit a quasi-static nucleation phase, observable quasi-static nucleation is not a prerequisite for dynamic rupture. Given an appropriate overstress, a sufficient impulse (provided here by a foreshock) can merge into dynamic rupture without an accelerating quasi-static slip phase. Similarly, although in our laboratory experiments all dynamic events exhibiting observable nucleation phases are initiated by seismic foreshocks, the impulse that initiates quasi-static nucleation can be supplied by any phenomenon causing localized acceleration, such as tectonic loading or other aseismic processes 20 . (2) How does the occurrence of a foreshock influence the nucleation of a large earthquake? On rate-weakening faults, foreshocks initiate a positive feedback loop between the sliding velocity and the stress drop, causing K bg to eventually grow to the dynamic range of the rate-dependent fault toughness K c (ν r ) with ν r ≈ c s . This drives the growth of the aseismically slipping patch and causes main-shock nucleation in sufficiently stressed conditions, making rate-weakening interfaces sensitive to foreshocks, especially when compared with rate-independent frameworks, and providing a mechanism compatible with a rate-dependent cascade up model 14 , even in cases in which foreshocks are the dominant driving mechanism of the rupture tip (that is, cascade-like). Larger foreshocks shrink the timescale and reduce the stress criticality required for nucleation. The largest foreshocks can themselves merge into a mainshock under sufficiently stressed conditions. Yet, observations of quasi-static nucleation phases in foreshock to main-shock sequences are relatively rare in nature. This is probably because the impulses of most foreshocks are too small to lead to nucleation after the initial slowdown (regime 1), too large to not directly merge into a mainshock (regime 3) or because quasi-static phases are simply too difficult to identify in observations 24 . (3) Does the characteristic slip distance representative of dynamic rupture also describe nucleation? These results allow for the inference of reasonable ranges (between 0.3 and 3.0 mm) of state-evolution slip distances, L, of real earthquakes during their nucleation phases (Fig. 4d). These are a notable departure from the larger values of L inferred for dynamic rupture; for example, the state-evolution distance of 0.8 m used for modelling Tohoku-Oki 23 .

Taken together, our results demonstrate that nucleation is a dynamic, rate-dependent process that can manifest in several forms and whose characteristics hold important implications for understanding the conditions under which rupture initiates and evolves. Beyond earthquakes, our results highlight a fundamental mechanism of frictional failure that may also govern the stability of engineered and tribological interfaces, as well as landslides and icequakes.

Methods

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 3 min −1 , 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,

(2)

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

(3)

in which ε v is the vertical direction and ε 315 and ε 45 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,

(4)

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 −1 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. 20 ), 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 fracDelta Tmu ^prime L, in which ΔT ≈ μ′δ a , such that the scaling parameter reads fracdelta _aL. 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. 46 )), 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. 14 , 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 20 ,

(5)

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 26 .

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 47 . Both methods are subject to several assumptions and potential errors. Illapel 2015 exhibited repeater earthquakes beginning approximately 140 days before the mainshock 9 . 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 4 . 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 3 , 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 40 . 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 10 . 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 2 . 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 48 . 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 6 . 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 13 . 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 49 ,

(6)

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 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 50 . 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 22 , 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 0 and ‘overstress’, the difference between the initial stress (ratio) f 0 = τ 0 /σ and the steady-state level of friction at the ambient sliding velocity, Δf 0 = f 0 − f ss (V 0 ). The fault is embedded in an infinite elastic medium and its rate-and-state-dependent resistance to sliding 32 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 16 or slip 52 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=-fracdtau ddelta), 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 _c=(fracDelta f_p2b)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 54 , 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 55 , 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=overlineK 2 /(2barmu ) into the propagating crack tip with the fracture energy of the frictional breakdown process G c . Here overlineK=k(nu _r)K is the dynamic SIF obtained from the static value K multiplied by the wave-mediated dynamic prefactor k(ν r ) and overlinemu =g(nu _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 22 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 22 gives

(7)

in which Li 2 is the dilog function. This expression further simplifies to G_c~ (1/2)(fracDelta f_pb) 2 bsigma 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. 22 . Here we use an explicit expression,

(8)

in which barnu 0 =exp (-fracDelta f 0 b)fracmu ^prime bsigma V 0 is the characteristic rupture velocity embodying initial fault conditions (ambient sliding velocity V 0 and initial overstress Δf 0 = f 0 − f ss (V 0 )), κ 0 ≈ 0.838 and ProductLog −1 (x) denotes the −1 branch of the Lambert productlog function providing the real-valued solution of ne n = 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 22 observes that Δf p (ν r ) is essentially indistinguishable for the two laws when plotted in his Fig. 4a.

Reasonable (zeroth-order) approximations, Delta f_p~ bln(nu _r/barnu 0 ) and G_c=(1/2)(ln(nu _r/barnu 0 )) 2 bsigma 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_c=(sqrtg/k)sqrt2mu ^prime G_c, with the prefactor (sqrtg/k)(nu _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,

(9)

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

(10)

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,

(11)

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

(12)

as parametrized here in terms of the fault initial state: f 0 , V 0 and overstress Δf 0 = f 0 − f ss (V 0 ). 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,

(13)

Alternatively, we rewrite

(14)

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. 22 ) and solving for V eff ultimately results in

(15)

Finally, the EoM

(16)

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_bg=Delta tau (V_eff)sqrtpi 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 0 defined as the distance between the background stress ratio f 0 from the steady-state friction value f ss (V 0 ) at the ambient sliding velocity V 0 —and transient rupture velocity ν r . In our experiments, the fault is probably neutrally or mildly understressed (Δf 0 ≱ 0; Fig. 3c and Extended Data Fig. 12). Furthermore, we introduce K_fs=Delta T/sqrtpi 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. 20 ) (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~ K_bg~ (b-a)sigma sqrtpi ell * ln(V_eff/V 0 ), balances the asymptotic expression for the fault toughness K_c~ sqrtmu ^prime bsigma L* ln(nu _r/barnu 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 sqrtpi ell ~ sqrtmu ^prime bsigma 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 0 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 0 and V 0 ) 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, fracDelta Tmu ^prime L, initial overstress, fracDelta f 0 b, ambient sliding velocity, fracV 0 (bsigma /mu ^prime )c_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. 56 ) while also performing a limited suite of solutions for the rock-like value a/b ≈ 0.75 (ref. 35 ) 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 fracDelta Tmu ^prime L (implicit), initial overstress fracDelta f 0 b and nondimensional ambient sliding velocity fracV 0 (bsigma /mu ^prime )c_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 < 0) faults when the foreshock force decreases to the Δf 0 -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, fracDelta Tmu ^prime L and fracDelta f 0 b, and ambient sliding velocity fracV 0 (bsigma /mu ^prime )c_s (Extended Data Fig. 9a). This collapse manifests two asymptotic behaviours in which Δt c ≈ L/V 0 is a constant for transients with very small minimum sliding velocity, comparable with the initial fault value, that is, when V eff,min ≈ V 0 , and Δt c ≈ (L/V eff,min )ln(V eff,min /V 0 ) when V eff,min ≫V 0 . Essentially, to the first order, Δt c is inverse in V eff,min with the low-velocity cut-off at about V 0 . Small departure from the collapsed curve on Extended Data Fig. 9a corresponds to the solutions with appreciable understress Δf 0 < 0. This second-order behaviour is more apparent when we remove the dominant (inverse in V eff,min ) behaviour when plotting the product Δt c V 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

(17)

with the transition function

(18)

between the small-velocity, mathcalD~ C 0 mathcalV, and large-velocity, mathcalD~ C 1 ln(fracmathcalVmathcalV 1 ), asymptotes, and fitting constants C 0 , C 1 , mathcalV 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 m in to the crack-average one V eff,min . To accomplish this, we use the EoM approximate solution for the sliding velocity distribution 22 V(x) = (V eff /2)(1 − (x/ℓ) 2 ) −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 56 , equation (17) was fit to the laboratory data, yielding a fault ambient sliding velocity V 0 = 1.35 × 10 −6 m s −1 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 56 .

Similarly, in Fig. 4d, equation (17) was used to prepare possible fits for the natural earthquake data assuming granite friction 35 with a = 0.014 and b = 0.019 and natural fault ambient sliding velocity V 0 = 10 − 12 m s −1 . 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 − 9 m s −1 . Thus, any suitable choice of the ambient sliding velocity assuming an initially ‘locked’ fault, that is, V 0 ≪ 10 − 9 m s −1 , 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 fracV eff,min V 0 ~ 10* e^fracDelta T0.51mu ^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 −1 , 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.

登录APP查看全文