Cascade and pre-slip models oversimplify the complexity of earthquake preparation in nature

Earthquake precursory processes have been central to scienti ﬁ c inquiry for nearly a century. Recent advancements in earthquake monitoring, geodesy, and data analysis including arti ﬁ cial intelligence, have substantially improved our understanding of how earthquake sequences unfold leading to the mainshock. We examine the available seismological and geodetic evidence describing preparatory processes in 33 earthquake sequences with M W [3.2 – 9.0] across different tectonic and stress conditions. Our analysis reveals common patterns, and sheds light on the interplay of structural, tectonic and other boundary conditions that in ﬂ uence the dynamics of earthquake sequences, and hence, in the seismo-geodetic observables prior to the mainshock. We place particular emphasis on connecting observed phenomena to the underlying physical processes driving the sequences. From our ﬁ ndings, we propose a conceptual framework viewing earthquake preparation as a process involving several juxtaposed driving physical mechanisms on different temporal and spatial scales, jointly leading to the stress increase in the future epicenter. A

A persistent goal in geosciences is to enhance protection against seismic hazard.To advance this topic, scientists are working to understand whether hazardous earthquakes exhibit some characteristic process, either at short or long spatiotemporal scale, prior to the mainshock [1][2][3][4] .Ideally, if such processes could be reliably and systematically observed prior to large earthquakes, automated near-real-time data workflows could enable extended warning times.However, detecting earthquake preparation processes, and understanding their significance remains a challenge in seismology.The current state of knowledge reveals a diverse set of observations, including sequences with no detectable preparatory processes.The tectonic conditions that facilitate or impede detectable preparatory processes before mainshocks are still not well understood 3,5 .Overall, various lines of research suggest that the prospect of short-term earthquake prediction, if achievable, remains a goal for the future of seismology 3,6 .
The fundamental physical and mechanical processes that lead to earthquake preparation and nucleation have been studied both theoretically and experimentally.In the laboratory, early work in the 1980's already achieved the observation of pre-failure processes 7,8 .These observations were supported by early fracture mechanics models 9 , which then evolved into widely accepted theoretical earthquake nucleation models, such as the slipweakening 9 and the rate-and-state 10,11 models.These models predict a quasistatic phase of deformation prior to earthquakes, during which the initially stable slip in the future plane begins to accelerate (pre-slip).This pre-slip continues until a critical nucleation size, characterized by a crack length L C , is reached 10,12 .At L C , the system becomes unstable and the rupture grows quasi-dynamically, i.e., an earthquake is initiated.A primary focus of these models and laboratory experiments has been to investigate the effect of stress and material heterogeneity within the pre-slip region in order to understand the occurrence of foreshocks 12 .These heterogeneities included structural fault properties such as roughness [13][14][15][16] , loading rate variation 17 , effective confining pressure 18 , normal stress 19,20 , fault strength 21 , temperature or the presence of fluids [22][23][24] .Another important focus has been to obtain direct observations of the microscale damage evolution during fault nucleation and propagation, and how they vary with respect to different properties, as for example rock type and stress 25 .A summary of the properties and processes that have been suggested to influence the preparation and nucleation of earthquakes is given in Box 1.
While laboratory earthquakes are, to some extent, understood, our comprehension of the processes leading to an earthquake in real faults is less clear [26][27][28][29][30][31] .To date, the upscaling of physical models and observations from the lab to nature remains a challenge.This is mainly due to the inherent complexity of major faults, and to our incomplete understanding of the physical processes that cause earthquakes in nature.In addition, the geophysical data available to characterize many earthquakes are limited and heterogeneous, especially when they occur in remote locations.As a result, many studies based on geophysical data show only sparse observations of the processes occurring prior to earthquakes, often leading to ambiguous interpretations of their preparatory processes (see e.g., the different mechanisms proposed for the 1999 M W 7.6 Izmit earthquake 29,32 , or discussion in Gomberg 5 . In this study, we undertake a comprehensive review of numerous studies that analyzed geophysical data to shed light on the collective patterns and physical processes influencing the occurrence of earthquakes in nature.Our analysis shows that (1) there is a multitude of processes that can lead to a stress increment on faults, eventually triggering an earthquake, and (2) these physical processes can occur concurrently at spatial and temporal scales Box 1 | Properties and processes influencing the preparation and nucleation of earthquakes

Fault zone material composition
Geodetic and seismological observations have documented or suggested the occurrence of slow slip and/or pre-slip preceding some mainshocks.The fault zone material composition strongly controls its frictional properties, which may promote (i.e., a velocity strengthening medium) or potentially inhibit (velocity weakening medium) slow slip 84 .For example, low-friction gouges rich in clay content have been observed to favor velocity strengthening 109,110 .In contrast, crystalline rocks were often found to be velocity weakening in the field, related to locked segments [111][112][113] .Note that velocity strengthening patches may also be able rupture dynamically 114 .In addition, a spectrum of rock damage characteristics depending on the deformed rock types were observed during the nucleation of cracks in laboratory faults using an X-ray triaxial apparatus 25 .These experiments suggested that changes in the local seismic velocity during nucleation are expected to vary according to the specific tectonic setting.For example, lower S-wave seismic velocities were resolved in a fault patch that hosted >20.000 foreshocks before a M 6 earthquake in the Gofar oceanic transform fault in the East Pacific Rise 115 .This fault patch acted as a barrier for the propagation of the mainshock, and the authors concluded that the rupture was controlled by the contrast of different rock materials along the fault.

Fault roughness
Fault roughness measures the deviation of a fault surface from planarity at all scales.Increased fault roughness is more representative of faults in nature and its effect on nucleation has been investigated in numerical simulations 15,40 and experimental work 14,[116][117][118][119] .In rock deformation experiments, slow slip and foreshock rates have been amplified by roughness 14 .In nature, however, direct measurements of fault roughness at seismogenic depths are typically not available.Recently, an approach that used the deviation of the seismicity from a planar fault surface to estimate the fault roughness has been proposed 120 .

Fault segmentation
Fault segmentation leads to more heterogeneous fault properties and stress field 121,122 .These structural heterogeneities strongly control the seismic energy release during the interseismic and coseismic periods.Increased geometrical complexities can also lead to the occurrence of slow slip 123 .In laboratory experiments, the occurrence of foreshocks and their spatial distribution is strongly controlled by the distribution of structural fault heterogeneities 12 .In nature, fault jogs, bends and step-overs have been observed to host numerous foreshocks for selected earthquakes in California 27,124 , in good agreement with the proposed models from the laboratory.

Pressure and temperature conditions
Increased normal stress (σ n ) suppresses small-scale fault heterogeneities, thereby inhibiting the occurrence of foreshocks 26 .Temperature is thermodynamically related to pressure.Early earthquake nucleation models 12 already suggested a strong variation in slip behavior as a function of pressure and temperature.The critical nucleation length L C has been shown to be inversely proportional to σ n in the slip weakening framework 125 .This inverse dependence between σ n and L C has been observed in friction experiments on polycarbonate surfaces 19 .Numerical simulations of faults described with a rate-and-state framework show that mainshocks tend to nucleate in the regions with higher σ n , representing asperities 15 .
Local σ n spectrum couples with fault roughness, generating an even broader of slip behaviors.Stick-slip experiments showed longer duration of slip events and lower stress drops at lower σ n than their counterparts at higher σ n .

Loading rates
Laboratory friction experiments on polycarbonates have demonstrated that L C is inversely proportional to the rate at which the system accumulates displacement (i.e., the loading rates) 17 .When the system is loaded more rapidly, it is more susceptible to dynamic rupture.The loading rate in the laboratory roughly corresponds to the tectonic loading rate in nature.In this respect, differences between foreshock sequences of M > 6.5 interplate earthquakes and their intraplate counterparts have been observed.Earthquakes from interplate regions where loading rates are typically larger exhibited accelerating seismic activity in the months to days prior to the mainshock, whereas this was not the case for intraplate settings 2 .The effect of loading rates has also been demonstrated in the behavior of repeater sequences during afterslip transients.Following the 2004 Parkfield M 6 earthquake, local earthquake repeating sequences displayed greatly reduced recurrence times that gradually increased thereafter 127 .Similar behavior has been observed in the creeping segment of the Main Marmara Fault 128 .
Presence of fluids and pore fluid pressure The presence of fluids can stabilize or destabilize faults.According to Terzaghi's principle, increasing pore fluid pressure reduces the effective normal stress, hence bringing a fault closer to failure and facilitating slip.Laboratory experiments on various rock types have shown that increased pore fluid pressure tends to promote velocity weakening frictional behavior rather than velocity strengthening, resulting in a higher probability of earthquake rupture 22,23 .In nature, fluid pressurization and migration along faults have been proposed to explain different types of seismic evolution.For example, the 2009 L'Aquila earthquake 89 , the 2016-2017 earthquakes in central Italy 92 , and the Iwate-Miyagi earthquake in Japan 129 have been suggested to occur in fault patches that were overpressured by fluids.
However, an increase in pore fluid pressure does not necessarily lead to dynamic rupture and has also been observed to play a key role in promoting slow slip on faults 130 .This behavior has typically been explained by a dependence of Lc on fluid pressure, which is explicitly assumed in both slip weakening and rate-and-state theories 131,132 .Recent experimental work has also shown that variations in the pore fluid pressure, which affect the initial effective normal stress, can generate the entire spectrum of slip behavior 133 .

Models describing earthquake preparation and nucleation
A combination of laboratory, field and theoretical studies led to the formulation of two end-member models of earthquake nucleation, which predict the spatial and temporal evolution of stress and strength in the vicinity of a mainshock.These models have been widely used to describe the geophysical observations prior to an earthquake and are summarized as follows.
(1) The pre-slip model predicts aseismic slip acceleration prior to an earthquake.The aseismic acceleration takes place within the nucleation zone, which increases in size until it reaches a critical nucleation length (L C ), from which the earthquake dynamic rupture begins 12 (Fig. 1a).It can be seen as a self-nucleation model, where the eventual observation of pre-slip indicates in advance that an earthquake will occur.Therefore, this model is compatible with earthquake predictability 33,34 .
One outcome of this model is a scaling between the space occupied by foreshocks and the size of the incipient mainshock 27,35,36 .Precursory aseismic acceleration has been observed in the laboratory 13 and in numerical models 12 .Experimentally, results consistent with this model have been observed in frictional rock deformation experiments on m-scale granite blocks with smooth faults performed at very low normal stress 20 .If present, aseismic pre-slip is expected to occur within the nucleation region of the future dynamic rupture, with the slip release accelerating as dynamic rupture approaches.In nature, however, there are almost no observations of earthquake mainshocks preceded by well-documented aseismic or slow deformation that unambiguously represents pre-slip.Earthquake sequences from tectonic environments such as oceanic transform faults, where aseismic deformation accommodates a substantial proportion of the plate movement, have been proposed to be particularly consistent with the pre-slip model 28,34 .
(2) The cascade model, in which the stress change associated with the occurrence of each foreshock triggers the next one, and eventually triggers the rupture of a larger earthquake by stress transfer (Fig. 1b).In this model, earthquake nucleation is viewed as a stochastic process in which each foreshock is no different from the mainshock.Therefore, the time and magnitude of the mainshock cannot be inferred before the end of the mainshock rupture 27 , giving no chance for earthquake predictability 37,38 .The physical models described above have often been proposed to explain the seismic and geodetic observations prior to earthquakes in nature 29,32,39 .However, observations from several sequences also suggest that these two models probably oversimplify the real processes.For example, pre-slip is mostly observed in homogeneous laboratory faults 20 , which are rarely found in nature.To account for the heterogeneous nature of real faults, new numerical studies 15,40,41 and laboratory settings 20,42 have been developed.The introduction of stress and strength heterogeneity in the models and experiments began to reveal complicated spatio-temporal patterns of foreshocks and aseismic slip, often resulting from stress feedback processes 15 , leading to new models describing these more complex deformation patterns.
(3) In the rate-dependent cascade-up model 20 , the nucleation is characterized by the coexistence of stress transfer (as in the cascade model) and aseismic slip, although the latter is not an intrinsic part of the nucleation process as in the pre-slip model.The heterogeneity in stress or strength strongly controls the nature of the observable processes before the mainshock, and the nucleation of an event is favored by an increase of stress in a rupture-prone fault region.( 4) The progressive localization model 43

Results: synthesis of observable processes preceding earthquakes
We compiled available studies of 33 earthquake mainshocks spanning a magnitude range M W [3.7-9.0] with the aim of summarizing the inferred physical processes that may drive their occurrence (Fig. 2a).The mainshock inventory includes earthquakes that occurred under various stress, kinematics, fault loading and structural conditions.Where possible, quantitative parameters describing the sequences were extracted from the reviewed articles (Supplementary Table S1).An individual description of each earthquake sequence is provided in Supplementary Note S1.Earthquake preparatory processes are mostly derived from seismic and geodetic observations, which are described separately below.

Seismological observables
In most of the studies reviewed, the assessment of mainshock preparation and nucleation was conducted by examining the spatial and temporal evolution of foreshocks.However, rarely were signals quantitatively linked to intrinsic nucleation processes.A possible exception is the M2016 M W 3.7 into Flats earthquake (Event #19 in Supplementary Note 1), where the mainshock was preceded by approximately 20 s of high-frequency radiation and a Very Low Frequency (VLF) earthquake, both located near the mainshock hypocenter 44 .
When comparing and reviewing different earthquake sequences, an important goal was to explore how the onset of an earthquake preparatory phase was defined.In laboratory experiments, this is typically associated with periods of visible inelastic behavior in the stress-strain curve 14,16 .Our study revealed that foreshock activity in nature was typically identified as anomalous with respect to the stable background seismicity, but identifying the onset of the stress-strain non-elastic behavior was rarely possible.Additionally, as source parameters (e.g., static stress drop, radiated energy, apparent stress) and the mechanics of foreshocks (e.g., fault orientation, frictional parameters) are a priori no different from mainshocks and aftershocks, there is no unique and unambiguous recipe for recognizing foreshocks before the mainshock occurs.For this reason, many studies defined the foreshock sequences as all the seismicity of smaller magnitude than the mainshock within a predefined time and space window 45 .The temporal scale typically ranged from minutes to months and the spatial scale  Circles below the panels represent events for which no foreshocks were identified and reported during the analyzed time period.A and B letters after specific earthquakes denote the corresponding catalog if more than one dataset was available (see Table S1 for details and references).
from meters to tens of kilometers around the mainshock hypocenter.The spatial and temporal windows selected in some of the early studies were probably designed to capture the short-term slip acceleration predicted by theoretical models 12 .Additionally, these choices were likely constrained by the absence of continuous geophysical time series covering annual periods until recent years.The application of more objective data analysis 43,46 , as for example, to the events preceding the 2014 Iquique earthquake 47 (see Event #24 in Supplementary Note S1) has helped to achieve a more rigorous estimation of the foreshock sequence.
In the following sections, we first analyze the influence of various monitoring setups and data analysis techniques on the seismic observables.Subsequently, we elucidate the connection between these seismic observables before mainshocks and the underlying physical processes.
Influence of monitoring conditions and data analysis.Ensuring that the monitoring conditions are sufficient to detect foreshocks is essential to correctly infer the physical processes driving the dynamics of the sequences.Our mainshock compilation covers the period from 1986 to 2023.During these years, the number of seismic stations available increased in many of the regions analyzed, resulting in an overall decrease in the magnitude detection threshold (see e.g., Hauksson and Shearer 48 ).In addition, development of advanced data analysis methods over the last decade, such as template matching 49,50 and machine learning 51 enabled further reductions in the magnitude detection threshold, typically by up to a unit of magnitude.The catalogs for the sequences reviewed were produced employing different data analysis techniques, and utilizing data recorded on seismic networks with different detection thresholds, hence the magnitude of completeness varies considerably (Fig. 2b, c).Differences in data quality between catalogs make it difficult to compare foreshock sequences objectively and may lead to differences in the foreshock rates.Furthermore, a high magnitude of completeness may obscure spatio-temporal foreshock patterns.Nevertheless, most monitoring setups capable of detecting seismicity at least three magnitude units lower than the mainshock detected foreshock activity 3 .
The analyzed time window preceding the mainshocks from the selected sequences spans from 30 min (for the 2019 M W 6.4 Ridgecrest 52 , Event #9 in Supplementary Note 1) to nine years (for the 2023 M W 7.8 Kahramanmaraş 53,54 , Event #3 in Supplementary Table 1).In Fig. 2b, c, we report the number of seismic events preceding the mainshock, as a function of the duration of the analyzed time period and the corresponding magnitude of completeness M C (where available), for seismicity catalogs derived using standard and enhanced techniques, respectively.The median M C of all the compiled catalogs derived by standard techniques is 1.5 (Fig. 2b), but for at least seven sequences (mainly before the year 2000) no M C was reported.In contrast, the median M C of the enhanced seismicity catalogs was reduced to 0.5 (Fig. 2c).The scatter in the observations is larger for seismicity catalogs derived by standard techniques (Fig. 2b) than those derived by enhanced techniques (Fig. 2c), revealing that standard data analysis is greatly affected by the local monitoring conditions.Despite a greater consistency between foreshock rates and analyzed time periods among the enhanced catalogs, there is still significant scatter (Fig. 2c).For example, looking at approximately one day before the mainshock, the number of reported foreshocks included in the enhanced catalogs vary from 50 for the 1999 Hector Mine (Event #7 in Supplementary Note 1) to almost 200 for the 2016 Te Araroa (Event #27 in Supplementary Note 1, Fig. 2b, c).
The variability of the observations presented here reflects to some degree the differences in station coverage and the methodologies employed to build catalogs.Nevertheless, most of the studies also presented quantitative measures of the temporal and spatial evolution of the seismicity (Table S1), associated with the physical processes which are described in Supplementary Note S1 and discussed in the following sections.
Foreshock dynamics from high-resolution catalogs.Upon minimizing differences in the monitoring conditions and data analysis among sequences, various seismic observables such as foreshock rates, moment release or the size of the foreshock area (Figs 2, 3) can be utilized to infer the physical processes in the lead up to an earthquake.In this regard, monitoring the time-space foreshock evolution is a common feature of the studied sequences 27,55 .Most of the studies analyzed reported that the evolution of foreshocks is connected to the complexity of the physical processes that preceded a mainshock.If foreshocks were a by-product of a single process, a specific time or seismic moment evolution would be expected, such as accelerated seismic release 37,56 similar to that observed in laboratory experiments 57 .However, a simple time evolution (e.g., acceleration) has rarely been observed.
The spatio-temporal migration of seismicity can be attributed to various physical mechanisms, such as aseismic slip 58 , fluid migration 59 or stress transfer between earthquakes 32 .Distinguishing between these mechanisms is aided by assessing the speed of migration 60,61 .Our analysis of foreshocks revealed a predominant intermittent temporal pattern (refer to Fig. 3a, b) accompanied by distinct spatial migrations (see Fig. 3c, d).These spatiotemporal patterns are interpreted as seismic manifestations arising from a diverse range of physical processes.For instance, foreshock sequences linked to the 2009 L'Aquila event (Event #1 in Supplementary Note 1) and the 2019 Marmara event (Event #31 in Supplementary Note 1) exhibit numerous temporal increments in seismicity (refer to Fig. 3a, b), with migrations partially indicating fluid propagation or aseismic slip, partly afterslip resulting from significant foreshocks, and partly indicative of stress interaction between asperities 62,63 .Furthermore, for nine of the sequences, foreshock spatiotemporal migrations were reported, with velocities ranging from 2 m/day to 20 km/day (see Supplementary Table 1).Migration velocities below 1 km/day are typically associated with fluid movement, as observed e.g., in the 2016 Pawnee (Event #14 in Supplementary Note S1, see ref. 64 ) or the 2010 El Mayor-Cucapah (Event #6 in Supplementary Note S1, see ref. 30 ).On the other hand, faster migration velocities are often linked to slow slip propagation, as exemplified in events like the 2016 Kumamoto (Event #8 in Supplementary Note S1, see ref. 55 ) or the 2014 Iquique earthquake (Event #24 in Supplementary Note S1, see ref. 65 ).
In summary, foreshocks exhibit an intermittent and spatio-temporally localized pattern, indicating an interplay of processes rather than directly tracking a single intrinsic mechanism leading to the mainshock.This ensemble of processes occurs over a wide area, as depicted in Fig. 3.A comparable complex behavior is commonly observed both experimentally and in numerical models describing earthquake ruptures on heterogenous faults.In experimental settings, the evolution of Acoustic Emissions (AE) recorded in rock deformation experiments before stick-slip events often showcases clusters of AE events localized in multiple asperities.These clusters tend to coalesce as the fracture approaches, reflecting a behavior similar to that observed in earthquake foreshocks 14,16 .
Can we learn about the upcoming mainshock size from the foreshocks?.A related topic of crucial importance, with implications for hazard mitigation, is whether the earthquake preparatory and nucleation processes can inform about the upcoming mainshock size.In the pre-slip nucleation model, foreshocks are a by-product of a growing slow slip, and they rupture small brittle asperities within the slow slip area (Fig. 1a).Accordingly, the area covered by the foreshocks could give an indication of the critical nucleation length L C of the mainshock, and hence, of its magnitude.Dodge et al. 27 estimated the area activated by the foreshock sequences of six strike-slip mainshocks (assuming a constant stress drop of 3 MPa) and compared the radius of the foreshock area with the mainshock coseismic moment.They found that the extent of the foreshock sequences, representing the size of the nucleation region, roughly scaled with the mainshock seismic moment.Chen and Shearer 60 further investigated this scaling for 14 foreshock sequences potentially driven by fluid, and concluded that there was no clear relationship.
We increased the number of observations with our additional sequences.Although this plot may be affected by the magnitude of completeness and the processing choices made by the different studies, it still shows that, the size of the area activated by the foreshocks does indeed tend https://doi.org/10.1038/s43247-024-01285-y to increase with mainshock size (Fig. 4), suggesting that aseismic slip may have been involved in mainshock nucleation at the spatial scale delimited by the foreshock sequence.However, several sequences for which an enhanced seismicity catalog is available do not follow such a trend.Among them, the foreshock sequences of strike-slip mainshocks such as the 1999 M W 7.1 Hector Mine (Event #7 in Supplementary Note S1), the 1992 M W 7.3 Landers (Event #5 in Supplementary Note S1) or the 2019 M W 6.4 first Ridgecrest (Event #9 in Supplementary Note S1), activated a smaller area than predicted by the overall trend (Fig. 4).These mainshocks invoke a potential rupture driven by more localized physical processes, such as static stress transfer (e.g., as discussed in Yoon et al. 36 for the 1999 M W 7.1 Hector Mine, Event #7 in Supplementary Note S1).Conversely, the foreshock sequences of the 2014 M W 8.1 Iquique (Event #24 in Supplementary Note S1) and 2017 M W 6.9 Valparaiso (Event #29 in Supplementary Note S1) subduction earthquakes activated a broader area than what predicted by this empirical trend (Fig. 4).This suggests that larger scale processes such as subduction slab-pull (see section "Larger-scale processes") also played an active role in controlling the evolution of the future earthquake region and prepared it for rupture.
Despite the potential trends between foreshock area and mainshock size mentioned above, the self-similarity of earthquakes in nature suggests that some ruptures that start as small earthquakes may grow larger, which is more consistent with the cascade model.Numerical simulations have shown that the final mainshock size is largely controlled by the dynamics of its rupture and propagation 66 .At the same time, numerical simulations and observations of highly correlated earthquakes in Japan have shown that the existence of hierarchical physics play an important role in controlling the rupture dynamics and hence, the final upcoming mainshock size 67,68 .

Geodetic observations of slow-slip preceding earthquakes
As with the seismological observations, the vast majority of geodetic studies focusing on the times leading up to large earthquakes showed a significant spatiotemporal complexity, including transient accelerations and afterslip induced by large foreshocks.If a single process with a specific temporal evolution was taking place prior to an earthquake (e.g., acceleration of slip 37 ), we would see its signature in the recorded geodetic deformation.However, none of the studies unambiguously reported an isolated acceleration behavior.Thus, they do not directly trace a process of intrinsic earthquake self-nucleation.
One of the most complex GNSS observations was reported before the 2011 M W 9.0 Tohoku-Oki earthquake (Event #21 in Supplementary Note 1).Prior to the mainshock, a slow transient was recorded in the GNSS instruments, beginning approximately 2 days before the mainshock and following a M W 7.3 foreshock 31,69 .The recorded aseismic deformation was initially interpreted to be aseismic pre-slip, and thus as evidence for selfnucleation 31 .Nevertheless, the pre-slip did not manifest within the anticipated future mainshock rupture zone.Additionally, the observed logarithmic temporal decay of the slip aligns with expectations for the afterslip following the M W 7.3 foreshock 70 .Inversion of the GPS data showed that this recorded aseismic slip preceding the mainshock was the most likely explanation for the occurrence of foreshocks and repeating earthquakes prior to the M W 9.0 31 .Therefore, this recorded aseismic slip (afterslip) probably does not represent pre-slip as an intrinsic nucleation process of the impending M W 9.0, mainly because it did not occur within the future mainshock rupture and decayed with time instead of accelerating.Nevertheless, this slow transient most likely contributed to the stress increase in the future mainshock area and thus to the nucleation of the 2011 M W 9.0 Tohoku-Oki earthquake 31,70 .
Another sequence highlighting the complex interplay between different slow slip processes prior to earthquakes emerged from the 2017 M W 6.9 Valparaiso earthquake [71][72][73] (Event #29 in Supplementary Note 1).The event was preceded by slow deformation recorded by geodetic instruments, which was further evidenced by the occurrence of repeating earthquakes 71,74 .This aseismic slip started after a M W 6 foreshock and decayed in time, again suggesting the occurrence of afterslip, as before the 2011 Tohoku-Oki earthquake.Caballero et al. 72 calculated the slip balance between foreshocks and aseismic slip, and concluded that seismicity and afterslip could account for only half of the observed slip.Moutote et al. 73 generated an enhanced catalog of seismicity framing the mainshock.They reported the occurrence of aseismic transients starting one day before the foreshocks and continuing for several days after the mainshock.The authors concluded that the seismicity from this sequence was driven by a slow slip transient that started before the mainshock and continued for several days after it.
Although most geodetic studies showed a complicated time evolution of slip prior to mainshocks, an intriguing new result comes from the joint   63 and Cabrera et al. 64 , with their inferred physical mechanism reported in the legend.
https://doi.org/10.1038/s43247-024-01285-yanalysis of GNSS data close to the future nucleation of 93 M > 7 earthquakes worldwide 75 .By stacking recordings from all these sequences, an accelerating signal starting 2 h before the mainshocks emerged, which could represent a potential precursory signal, intrinsically linked to the nucleation process.

Larger-scale processes
With the collection of decades of geophysical data over entire plate boundary regions, scientists began to look for larger-scale and longerterm signatures of processes leading to large earthquakes.From the analysis of seismic catalogs 76 , geodetic time series 77 and gravity anomalies 78 in Japan and Chile at the plate boundary scale, authors reported a long-term (years) slab-pull effect favoring the occurrence of large earthquakes.However, the subtle signals from these studies have been a matter of debate and reported to not be significant by other authors 79,80 .Other studies highlighted how stress transfer from distant earthquakes (100-1000 km) can affect the fault coupling in some regions of the megathrust and promote the weakening of the interface 71,74,81 .Such behavior was recently reproduced in laboratory analog experiments recreating subduction zones 82 .Earthquake preparation on large spatio-temporal scales has also been investigated for strike-slip faults.Ben-Zion and Zaliapin 43 utilized the decades-long Southern California seismicity catalog to show that four of the largest recent earthquakes in California, were preceded in the previous decades by generation of off-fault rock damage around the eventual rupture zones.Seismicity localization around the main fault segments starting 2-3 years before the mainshocks was also reported.

Summary: a conceptual multi-scale model for earthquake preparation and nucleation
From our review analysis, it emerges that the seismic and geodetic observations that precede mainshocks can be affected by the resolution of the monitoring setup and the data processing techniques.Furthermore, the choice of temporal and spatial windows of analysis imprints an upper boundary to the characteristic scales of the physical processes that are possible to resolve.Beyond these differences, we posit that these studies also reveal a genuine diversity of physical processes with different characteristic scales that increase the stress on a fault and ultimately bring it to failure, acting individually or juxtaposed.
Among the reviewed studies, almost no sequence reported observations that can be directly related to the earthquake self-nucleation (with possible exceptions being ref. 44on the Minto Flats earthquake, see Event 19# in Supplementary Note S1 and Bletery and Nocquet 75 in stacks of several sequences as discussed in Section "Geodetic observations of slow-slip preceding earthquakes").Most of the slow slip observations reported so far in different studies did not represent an intrinsic part of the nucleation.This suggests that, if self-nucleation processes as those resolved in laboratory experiments occur at larger scale in nature, their detection is very difficult 83 .This may be because their associated signals are likely to be too small in amplitude or too short in time, and are often convolved with other processes, which prevents them from being systematically identified by our current instrumentation and data analysis techniques.
Foreshocks are the most common observation preceding earthquakes 26,27 , but most studies support the notion that foreshocks are a by-product of processes that alter the state of stress in a given fault and reflect its heterogeneities.Only rarely we observe clear foreshock accelerations as predicted by numerical modeling or in some laboratory experiments 10,14,15,20,[84][85][86][87][88] .When it does occur, acceleration is typically observed on a regional scale rather than being part of a localized process 53 .More often, foreshocks occur in spatiotemporal clusters likely to represent creep surges 15 , fluid pressure fluctuations 89 , or temporal clustering associated with static and dynamic triggering 32,36 .In the laboratory, analogous complexities in the observations and feedbacks between different processes are more evident at fractured fault interfaces with many asperities of variable roughness 90 .
In most studies of the sequences analyzed, a shear stress increment in the future mainshock rupture area was inferred, either by Coulomb stress transfer analysis or Brune stress drop calculation or through monitoring of indirect stress proxies, e.g., by resolving a temporal decrease in the Gutenberg-Richter b-value.This stress increase can be due to a large variety of physical processes (some of the most important listed in Fig. 5), as for example stress transfer from previous seismicity 32 , the occurrence and propagation of slow slip transients 91 , the migration of fluids along a fault patch resulting in a pore fluid pressure increase 92,93 or any combination of these and additional mechanisms.Each of these mechanisms can operate at a different spatial scale, and they can also temporally vary.The spatial and temporal scales at which each of these processes may operate are influenced by the heterogeneities of the medium and connected to the fault geology, rheology, geometry, and segmentation.Larger heterogeneities in the medium are translated into corresponding heterogeneities in the state of stress at different scales.
Taken together, the observations from our review emphasize the distinctiveness of each mainshock, marked by its unique stress history across various scales, and leading to slightly different preparation phases for each event.Considered alongside variations in the resolution and magnitude of completeness in different studies, these disparities offer an explanation to the lack of scaling observed in several parameters extracted from the reviewed studies (Supplementary Table S1).Building on these findings, we propose a conceptual model in which the preparation and nucleation of earthquakes are shaped by a diverse array of physical processes, sometimes acting juxtaposed.These processes operate on different temporal and spatial scales, collectively contributing to stress build-up and fault slip (refer to Figs 5 and 6).
Experimental work, numerical modeling and field observations indicate that the earthquake preparatory process includes damage localization, coalesce of fault segments leading to the failure of larger asperities and enhanced interaction of seismic and aseismic deformation 4,15,25,43,94,95 .Each of these processes leads to variations in the stress field operating at different spatial scales, and reflecting different heterogeneities and temporal evolution over the fault zone.
A multi-scale model for earthquake preparation is also supported by experimental data.Acoustic emissions recorded during rock deformation experiments on rough faults documented that the evolution of the employed https://doi.org/10.1038/s43247-024-01285-yseismo-mechanical parameters (e.g., b-value, clustering properties, event proximity, fault plane and stress variability, etc) displayed different trends depending on the different spatio-temporal scales considered.This supports the notion of a diversity of physical processes acting juxtaposed at different scales during the earthquake preparatory process 95 .As stress heterogeneities over various lengths are observed to increase with increased roughness and fault heterogeneity, our model might be more representative for seismic sequences occurring on juvenile or segmented fault structures.
The evolution of stress at the future mainshock hypocenter can be additionally influenced by other large-scale external processes such as seasonal or multiannual stress variations 96 and tidal stress evolution [97][98][99] which can either promote or inhibit the stress build-up on the fault, and hence, trigger or delay the nucleation of earthquakes.Finally, anthropogenic activities such as subsurface mining, fluid injection and extraction, or reservoir impoundment can also perturb the stress field over a wide spatiotemporal scale from the source region.These additional sources of stress perturbation add to those discussed in our work, and further complicate the time-dependent stress state of a fault (Figs 5, 6).

Research perspectives
Decades of research into the nucleation and preparation of tectonic earthquakes suggest that a multitude of physical processes operating on different temporal and spatial scales can influence the stress state of a fault, ultimately leading to its failure (Fig. 6).The interplay of different physical processes, combined with heterogeneous local and regional monitoring conditions, control which seismic and geodetic observables precede mainshocks.Recent detailed studies of observable slip processes before earthquakes represent a significant advancement in documenting the complexity inherent to the development of each earthquake sequence.These studies can now serve as examples for a more systematic evaluation of future mainshock initiation.To gain new insights on the preparation and nucleation of earthquakes, a collaborative effort is essential, involving laboratory experimental studies, the development of theoretical models and the monitoring of observables in natural settings.Here we summarize some of the relevant aspects of this collective endeavor.
(i) Densification of near-fault monitoring.Densification of near-fault seismic and geodetic monitoring is the most pressing need to gain insights into the fault activity prior to mainshock ruptures 100 .Improving the resolution can enable the detection and characterization of smaller foreshock activity, or resolve subtle geodetic signatures of aseismic processes or fluid movements, thus shedding more light on the underlying mechanisms 3 .An optimal seismic monitoring network can greatly contribute to the elucidation of fault behavior prior to a large earthquake, such as in the 2023 M W 7.8 Kahramanmaraş earthquake 53,54 .As we are witnessing that large earthquakes nucleate not only along main fault strands but also on secondary branches, these efforts should not be limited to main fault zones.(ii) Next-generation observations of fault processes over different frequency bands.Along with the need for instrumental development, there is a need to incorporate cutting-edge data.So far, we have seen how enhanced seismicity catalogs with a lower magnitude of completeness effectively illuminate what happens before the mainshock rupture (i.e., Fig. 2).A necessary step is to continue to generate enhanced observations for past and future earthquakes spanning different magnitudes, faulting styles and background stress conditions, with the aim of extracting similarities and differences between sequences.To this end, data analysis techniques employing artificial intelligence have a tremendous potential to boost fault signal observations 51,101,102 , which, in turn, may lead to clearer observations of earthquake preparatory processes.(iii) Automated multi-parametric and multi-scale workflows.With highresolution data and cutting-edge methods there will be the need to develop and test automated workflows based on time series of parameters 103 or data features [104][105][106] covering a wide range of spatiotemporal scales.Such workflows could shed light on which physical processes are most important in the preparation and nucleation of each mainshock, and at which scales these processes are most resolvable.They could also help to distinguish between the different juxtaposed physical processes operating at different spatiotemporal scales.(iv) Identifying the onset of the preparatory phase.Experimentally, the onset of the preparatory phase is typically considered to be the time when the stress-strain curve leaves the elastic regime and the specimen deforms non-linearly.A promising research target is to determine whether the observables that we are able to monitor on faults in nature (e.g., clustering of seismicity, geodetic signals, increased energy release) are anomalous with respect to the stable conditions, and thus, represent part of the preparatory stage.For example, to be able to distinguish whether a burst of earthquakes is unlikely to continue (forming a socalled swarm) or, on the contrary, it could destabilize a fault towards a large rupture.To this end, underground laboratory experiments under controlled conditions have the potential to bridge the gap between the controlled conditions of laboratory experiments and the greater heterogeneity and complexity of faults in nature.For real faults, catching the evolution of susceptibility of earthquake triggering to periodic stresses (e.g., tides) might help to better identify those times in https://doi.org/10.1038/s43247-024-01285-ywhich a fault enters a critical state and is getting closer to the rupture, as recently suggested by laboratory 107,108 and field works 98,99 .(v) Revised physical frameworks describing earthquake nucleation.Some well-documented processes preceding failure in experiments and nature suggest an interaction between localized seismicity clusters around the future mainshock, and a complex interplay between seismic and aseismic deformation.This behavior is not yet well reflected in commonly used physical frameworks describing earthquake

Review article
nucleation such as the rate-and-state or slip weakening models.Therefore, there is also a need to revise these frameworks according to the new observations, perhaps allowing for a greater level of complexity or larger spatial and temporal scales.
emerged from recent geophysical observations over longer spatial and temporal scales before large earthquakes, revealing a somewhat larger variability than that captured by the theoretical and laboratory models.The progressive localization model describes a gradual evolution from distributed damage in a rock volume to more localized slip.During this localization process, many clusters of seismicity are expected to be observed within a zone containing multiple faults at different scales.Each cluster may have its own foreshocks, one of which may lead to the initiation of the mainshock rupture.In this model, the presence of an existing fault that is weaker than the surrounding medium is less (or not) relevant43 .(5) With the aim of generalizing the various existing models, an integrated model 4 was proposed describing the generation of large earthquakes.The initial phase corresponds to the localization model, with progressive generation of rock damage and ongoing seismicity along the future rupture zone.Foreshocks are driven by the stress transfer from the occurrence of previous seismicity or the presence of aseismic slip.Then, as in the rate-dependent cascade-up model, a foreshock can trigger the large dynamic rupture.

Fig. 1 |
Fig. 1 | The classical models for earthquake nucleation (after Ellsworth and Beroza 35 ). a Pre-slip model, (b) Cascade model.L C : critical nucleation length.The green outline represents the slip front.

Fig. 2 |
Fig. 2 | Earthquake locations and number of events from the compiled datasets.a Location of the earthquake mainshocks (stars) included in this analysis.Color is encoded with faulting style, with red, green and blue denoting normal, strike-slip and reverse faulting, respectively.Number of seismic events preceding the mainshock versus analyzed time span using catalogs derived with standard (b) and enhanced (c) data analysis techniques, respectively.Colored circles and squares denote sequences where M C or M min was available, respectively, where the color is encoded with M C or M min .White circles are otherwise utilized.Symbol size is encoded with magnitude.Circles below the panels represent events for which no foreshocks were identified and reported during the analyzed time period.A and B letters after specific earthquakes denote the corresponding catalog if more than one dataset was available (see TableS1for details and references).

Fig. 3 |
Fig. 3 | Foreshock evolution as a function of time and distance to mainshock for the 2019 M W 5.8 Marmara (Event #31 in Supplementary Note S1) and 2009 M W 6.1 L'Aquila sequences (Event #1 in Supplementary Note S1).a, b Inter-event time as a function of time to mainshock (gray circles) and cumulative number of events (purple line).Concentrations of gray circles reveal seismicity accelerations.c, d Distance as a function of time to mainshock.Spatial clustering is also observed during periods of accelerated seismicity.Colored boxes mark seismicity clusters discussed in Durand et al.63 and Cabrera et al.64 , with their inferred physical mechanism reported in the legend.

Fig. 5 |
Fig.5| Physical processes at different spatial and temporal scales.Chart summarizing some of the identified physical processes (in black Roman) and seismological and geodetic observables (in blue italics) that occur prior to mainshocks at various spatial and temporal scales.

Fig. 6 |
Fig. 6 | Conceptual sketch of physical processes influencing earthquake nucleation.Conceptual sketch illustrating some of the physical mechanisms at different spatio-temporal scales that may be involved in the preparation and nucleation of an earthquake.Temporal scales represent approximately (a) from hours to days, (b) from days to months and (c) from months to years.
36timated radius of the epicentral foreshock area versus mainshock coseismic moment.Color is encoded with the mainshock kinematics.Solid and dashed orange lines represent a linear fit to the data and one standard deviation, respectively (after Dodge et al.27; Chen and Shearer 60 ; Yoon et al.36).