Earthquake Nucleation Along Faults With Heterogeneous Weakening Rate

Abstract The transition from quasistatic slip growth to dynamic rupture propagation constitutes one possible scenario to describe earthquake nucleation. If this transition is rather well understood for homogeneous faults, how the friction properties of multiscale asperities may influence the overall stability of seismogenic faults remains largely unclear. Combining classical nucleation theory and concepts borrowed from condensed matter physics, we propose a comprehensive analytical framework that predicts the influence of heterogeneities of weakening rate on the nucleation length Lc for linearly slip‐dependent friction laws. Model predictions are compared to nucleation lengths measured from 2D dynamic simulations of earthquake nucleation along heterogeneous faults. Our results show that the interplay between frictional properties and the asperity size gives birth to three instability regimes (local, extremal, and homogenized), each related to different nucleation scenarios, and that the influence of heterogeneities at a scale far lower than the nucleation length can be averaged.

2 of 11 with depth due to changes in the local host rock lithology, roughness, or in situ conditions (normal stress, temperature, pore fluid pressure, etc.) (Ohnaka, 2003;Tse & Rice, 1986). Then, how do these multiscale heterogeneous frictional asperities influence the global stability of seismogenic faults? Recent studies (Albertini et al., 2020;de Geus et al., 2019;Dublanchet, 2018;Perfettini et al., 2003;Ray & Viesca, 2017, 2019 provide valuable insights on how heterogeneities impact the overall stability of frictional interfaces, but arguably oversimplify the complexity of natural faults by assuming either a homogeneous weakening rate (Albertini et al., 2020) or orderly placed asperities of uniform size (Dublanchet, 2018;Perfettini et al., 2003;Ray & Viesca, 2017, 2019. A comprehensive framework, which links the variations of frictional properties at all scales to the overall fault stability, is thus dearly lacking. In this Letter, we build on the theory of static friction (Rubin & Ampuero, 2005;Uenishi & Rice, 2003;Viesca, 2016a) and the physics of depinning (Cao et al., 2018;Démery et al., 2014;Tanguy & Vettorel, 2004) to develop a theoretical framework that predicts, for any heterogeneous linearly slip-dependent fault interface, the critical size c of the earthquake nucleus. Supported by numerical full-field dynamic calculations, we show that the nucleation of an earthquake is not only always triggered by the weakest heterogeneity, but can also emerge from the collective depinning of multiple asperities. We highlight that this shift in instability regime stems from the interplay between the characteristic size of the heterogeneity and the length scale set by the distribution of frictional properties. Finally, we show that, in assessing the stability of an interface, one has to mainly account for perturbations whose wavelength exceeds the nucleation length, since the influence of small-scale asperities can be averaged.

Dynamic Simulations of Earthquake Nucleation Along Heterogeneous Faults
We consider two homogeneous 2D semi-infinite elastic bodies that are kept in contact with a uniform normal pressure n , idealizing the fault structure as a planar 1D frictional interface indexed by . The fault is loaded through a macroscopic shear stress ∞( ) that slowly increases in time . The friction f that opposes interface motion is assumed to be linearly slip dependent and fluctuates along the fault (Figure 1a). It locally evolves as slip grows from its peak value p( ) to its residual one r ( ) with a weakening rate ( ) that describes the material brittleness/ductility. Variations of the frictional properties ( p, r , ) may emerge along natural faults due to local changes in geometry, roughness, lithology, or ambient conditions (Ohnaka, 2003;Tse & Rice, 1986). Recent works show that the nucleation along homogeneous (Viesca, 2016a) and heterogeneous (Ray & Viesca, 2017) faults in the (aging) rate-and-state framework could be investigated from the stability of an equivalent interface with spatially dependent piecewise linear slip-weakening friction. Despite restrictive assumptions, our work may then provide ways to predict rupture nucleation for more complex and experimentally supported friction laws.
As the macroscopic loading ∞( ) grows, it locally exceeds the friction p( ) , and the two bodies detach from one another by a slip ( ) (Figure 1b). Provided that the fault has been at rest for a time far larger than that set by the propagation of elastic waves, the evolution of is described by the quasi-dynamic equations of elasticity for Mode II cracks (Lapusta et al., 2000;Rice, 1993): where ∞ is the far-field macroscopic loading, s the shear wave velocity, * = ∕(1 − ) ( and being the shear modulus and the Poisson's ratio, respectively), and ′ is a linear operator.
In Equation 1, the term − * 2 s , often called "radiation damping," physically represents wave radiation from the interface to the two elastic bodies, while * [ ] represents the nonlocal contributions of the overall slip to the local stress state. To investigate the stability of such a heterogeneous fault, we run periodic dynamic simulations building on a spectral boundary integral formulation of fracture (Breitenfeld & Geubelle, 1998;Geubelle & Rice, 1995). These simulations account for both the static redistribution of stress of Equation 1 and dynamic stress transfers (see Section 1.2 in Supporting Information S1).
How is the fault stability influenced by the steadily increasing loading? It results in rather complex dynamics as can be observed in Figure 1b. Multiple regions slipping at an accelerated rate, referred to as "slip patches," start to nucleate on the positions where p is low. As the loading is further increased, they grow quasistatically and coalesce into larger slipping regions. This initial nucleation stage of duration Δ nuc proves rather quiescent since no major velocity burst is observed. Yet, at = 0 , an instability develops on the right part the fault: a rupture propagates dynamically, and the two bodies start sliding one onto another at a uniform slip rate.
If such a simulation constitutes one realistic scenario for natural earthquakes nucleation, the simultaneous growth of multiple slip patches prevents any accurate measurement of the size c of the instability nucleus, which might well be twice as large as our measurement of Figure 1c. Yet, identifying this critical length scale proves crucial since it gives access to (a) the loading levels (Uenishi & Rice, 2003) and (b) the position at which an earthquake nucleates (Albertini et al., 2020;Ampuero et al., 2006) when is homogeneous along the fault.
These difficulties arise from spatial variations of peak strength that have been proven to play no role in the stability behavior of a slip-dependent frictional interface, which is solely controlled by the weakening rate (Favreau et al., 1999;Uenishi & Rice, 2003). Indeed, assuming that the macroscopic loading ∞ slowly increases enough and that the slip perturbation is small enough, the interface velocity =̇ is described near the instability by Uenishi and Rice 2003: where only is involved. This observation is supported by recent numerical simulations of crack nucleation along interfaces with stochastic distributions of p and homogeneous (Albertini et al., 2020), except in rare situations where the asperity scale interacts with the nucleation length (Schär et al., 2021). One may Figure 1. (a) Nucleation dynamics along a 1D heterogeneous coplanar fault constituted of brittle (in orange) and ductile (in black) asperities; inset: the fault frictional properties locally follow a linear slip-dependent law; the frictional stress f of the interface goes from its peak value p to the residual one r when the local slip reaches its critical value c , defining a weakening rate = ( p − r )∕ c . p , r , and are varying independently along the fault position . (b) The heterogeneous fault is subjected to a uniform shear loading ∞ . Under the influence of the steadily increasing loading, several regions of the fault start slipping where ∞ locally exceeds the strength p . (c) Slip growth develops quasistatically without any significant velocity burst, until one slip patch reaches a critical length c that leads to the dynamic rupture of the whole interface (see Movie S1).

Geophysical Research Letters
LEBIHAIN ET AL.

10.1029/2021GL094901
4 of 11 then focus on variations of weakening rate to quantify the influence of multiscale heterogeneities on fault stability.

Measuring the Nucleation Length in Presence of Weakening Rate Variations: A Model Fault Approach
We thus focus on the stability behavior of an idealized fault along which both the peak p and the residual friction r are uniform ( Figure 2a). To make any parallel to Mode I fracture easier and without any loss of generality, we set r ( ) = 0 (Albertini et al., 2020). Meanwhile, the weakening rate may vary from several orders of magnitude along the fault. Following the procedure of Albertini et al., 2020 (see Section 1.1 in Supporting Information S1), we generate fields that follow Gaussian correlations up to a characteristic length scale . Moreover, the values of follow a beta distribution of average ⟨ ⟩ and standard deviation w between two extremal values [ min, max] . We set the nucleation length of the reference homogeneous material with uniform ⟨ ⟩ as the adimensionalizing length of the system hom c ≃ 1.158 * ∕⟨ ⟩ (Uenishi & Rice, 2003). In the following, we consider w = ⟨ ⟩ , min = 0.25⟨ ⟩ , max = 4⟨ ⟩ , and = 0.05 hom c . The behavior of such a heterogeneous interface remains out of scope of the current theories of rupture nucleation where is homogeneous (Albertini et al., 2020;Ampuero et al., 2006;Favreau et al., 1999;Uenishi & Rice, 2003). We then wonder how local variations of as well as their intensity impact the overall fault stability. p and r are considered uniform along the interface, while the weakening rate varies with the position. These variations occur over a characteristic length and are distributed following a beta distribution p( ) between two extremal values [ min, max] . (b) The model interface is loaded through an overstressed patch that slowly expands in time. A slip perturbation and an associated velocity perturbation =̇ develop as time grows. The dynamics are characterized by two phases: the first phase consists of quasistatic growth (solid lines) and the second involves dynamic crack propagation (dashed lines) when the slipping region reaches a critical size c . (c) This shift in dynamics is observed on the temporal evolution of the slip patch size ( ) or its growth velocity ̇ . In particular, ̇ hits an inflection point at = c (in linear-log space), which provides an accurate measurement of c from the growth rate ̈∕̇2 . See Movie S2, Movie S3, and Figure S2 in the Supporting Information S1 for comparison with the homogeneous case of Uenishi and Rice (2003).

Geophysical Research Letters
In presence of spatial variations of , the nucleation length is expected to fluctuate along the fault. In order to investigate the local instability dynamics, we force nucleation at a given point, referred to as "fault center," by considering a macroscopic loading consisting of a slowly expanding region of size = ( ≪ s ), where the stress locally exceeds the frictional resistance ∞ ≥ p (Figure 2b). We observe in Figure 2b that a typical nucleation event is very similar, yet much simpler, to that of the heterogeneous fault of Section 2.1. Its dynamics consists of two distinct regimes: (a) the first regime involves a stable quasistatic slip growth for 0 where a portion of the interface is slipping, while (b) the second involves an unstable dynamic crack propagation for 0 where a rupture front propagates until the whole fault is moving. The shift from one regime to another occurs when the slipping region outgrowths a critical length c independently of the nature of the loading shape as long as it is peaked (see Figure S4 in Supporting Information S1). Importantly, this instability results from the collective motion of multiple asperities ( c ≃ 13 in Figure 2).
To quantify the influence of spatial variations of on c , we first propose a heuristic framework to measure it from numerical calculations. Looking at the evolution of the slip patch size over time in Figure 2c, we observe that its growth velocity ̇ follows an S-shaped curve and hits an inflexion point (in linear-log space) when = c , as previously observed in laboratory experiments of earthquake nucleation between two polycarbonates blocks (Latour et al., 2013). The nucleation length c may then be estimated from the maximal growth rate ̈∕̇2 , as the length where the patch expansion is at its strongest. The validity of our heuristic approach is assessed on homogeneous interfaces for which hom c ≃ 1.158 * ∕⟨ ⟩ is known a priori (Favreau et al., 1999;Uenishi & Rice, 2003). We use it to numerically estimate the critical length c of heterogeneous interfaces with a ±5% precision, corresponding to the error observed for the homogeneous interface of known c (see Section 2.1 in Supporting Information S1).

Theoretical Model for Nucleation Length Predictions
Here, we propose a way to determine the critical length c analytically building on both the theory of static friction and the physics of depinning. For a velocity perturbation centered in = 0 with a support of size , Equation 2 becomes: * ′ is the linear operator introduced in Dascalu et al. (2000) and Uenishi and Rice (2003), and = 2( − 0)∕ is the reduced position.
To assess the fault stability, we perform a Linear Stability Analysis on Equation 3. It consists of finding the perturbation size for which the linear symmetric operator [ ]( ) = 2 * 1[ ]( ) − ( ∕2) ( ) admits a zero eigenvalue. Condensed matter physics provides one way to tackle this problem in the presence of heterogeneities (Cao et al., 2018;Démery et al., 2014;Tanguy & Vettorel, 2004): We expand the perturbation with the disorder intensity w up to the second order = 0 + w 1 + 2 w 2 and solve the eigenproblem [ ] = , where = 0 + w 1 + 2 w 2 (see Section 3 in Supporting Information S1). The nucleation length c is the solution of the transcendental Equation 4, which encompasses the main novelty of the paper. 10.1029/2021GL094901 6 of 11 (Albertini et al., 2020;Uenishi & Rice, 2003). The third term corresponds to the second-order contributions up to a critical mode c ≃ 2 c∕ . This higher order term accounts for the influence of the spatial shape of in all its complexity, beyond the special cases of periodic ordered distributions of asperities (Dublanchet, 2018;Perfettini et al., 2003;Ray & Viesca, 2017, 2019. Equation 4 is only valid as long as no point of the fault reaches its residual friction value (i.e., ( ) < c( ) ). But ( can be replaced by the instantaneous weakening rate [ ( )] to assess fault stability around a stable slip state ( ) in the case of nonlinear or piecewise linear slip-dependent friction. Yet, if Equation 4 gives qualitative information on the influence of a nonstationary weakening onto the nucleation process, it does not provide a quantitative framework to predict rupture nucleation for these more complex friction laws, as the slip evolution remains unknown.
When we compare the theoretical predictions pred c of Equation 4 to the numerically estimated critical lengths meas c for an asperity size that varies over four orders of magnitude, we observe an excellent agreement (Figure 3a). Note that (a) the nucleation length cannot be estimated from Uenishi and Rice (2003)'s homogeneous theory (Figure 3b) and (b) the second-order contributions are required for accurate predictions (see Figure S7 in Supporting Information S1).
Equation 4 unveils rich physics about the impact of microscopic heterogeneities on the macroscopic fault stability. From it, one can directly link the spatial profiles of weakening rate to the local evolution of the nucleation length c along the interface (see black solid lines in Figures 3c-3e). In our simulations, the effective nucleation length corresponds to the one predicted at = 0 due to the peaked nature of the loading. In more realistic cases, the position 0 of the earthquake nucleus (and the associated nucleation length c( 0) ) will depend not only (a) on heterogeneities of peak strength p , but also (b) on the spatial shape of the macroscopic loading ∞ , similarly to what has been observed for nonlinear slip-weakening laws (Rice & Uenishi, 2010).
Overall, our framework provides ways to quantify the influence of a single heterogeneity on the fault stability depending on its size and intensity (see Figure 3) as well as that of the superposition of multiple perturbations of frictional properties (see Figure 4). Next, we build on Equation 4 and distinguish three instability regimes that can be linked to realistic earthquake nucleation scenarios on natural faults.

Instability Regimes in Earthquake Nucleation
In natural fault zones, heterogeneities in friction occur over many different scales. We observe them at the centimetric scale with minerals, clasts, and foliation, at the metric/decametric scale along large faults consisting of different lithologies, up to the scale of tectonic plates where kilometric asperities generated by a heterogeneous stress distribution have been suggested as potential nucleation sites for megathrust earthquakes in subduction zones. It is still uncertain how those different scales may interact with each other and how they ultimately impact the nucleation of earthquakes. Building on Equation 4, we highlight in Figures 3c-3e three different instability regimes, referred to as local, extremal, and homogenized regimes, respectively. They emerge from the interplay between three length scales: the heterogeneity size , the nucleation length associated to average frictional properties hom c ≃ 1.158 * ∕⟨ ⟩ , and the scale set by the weakest defect along the fault min c ≃ 1.158 * ∕ max .
1. Local regime: when ≳ hom c and min c , the weakening rate is almost constant over the nucleation patch (see Figure 3c). One then retrieves the dynamics of homogeneous nucleation (Uenishi & Rice, 2003), and the effective nucleation length is set by the local frictional properties at the fault center c(0) ≃ loc c ≃ 1.158 * ∕ (0) , which can be distributed above ( (0) < ⟨ ⟩ ) or below ( (0) > ⟨ ⟩ ) hom c 2. Extremal regime: when ≪ hom c and ≳ min c , a critical nucleation patch of size c(0) ≃ min c may develop within a single brittle asperity of size , where the weakening rate reaches its maximal value max . This small event destabilizes the interface as a whole, generating a complex dynamics of multiple slip pulses (see velocity map in Figure 3d for which = 0.01 hom c = 2 min c ). Along natural faults, these small ruptures may be arrested by local barriers of strength p , but they may trigger a cascade of nucleation events centered on other weakest spots until the entire fault fails (de Geus et al., 2019;Noda et al., 2013;Zhang et al., 2003). Note that (a) these brittle asperities influence the effective nucleation length c far

Geophysical Research Letters
LEBIHAIN ET AL.
10.1029/2021GL094901 7 of 11  Uenishi and Rice (2003) loc c ( ) ≃ 1.158 * ∕ ( ) at = 0 . The interplay between the length scales set by the frictional properties and the asperity size gives birth to three instability regimes: (c) When is larger than the homogeneous nucleation length

Influence of Each Asperity Scale to the Global Stability of Heterogeneous Fault
So far, we considered cases where the distribution of weakening rate asperities could be described through a unique length scale . Yet, the heterogeneities of weakening rate may emerge from, for example., the fault roughness that exhibits a scale-free self-affine behavior that spans over several decades of length scales (Candela et al., 2012), which makes the modeling of rough faults particularly challenging from a numerical point of view. Up to now, it is still largely unclear which length scales actively participate in the fault stability and which may be averaged in a realistic modeling of earthquake nucleation along rough faults. The influence of the superposition of a unimodal perturbation ̂ of period pert and amplitude pert to the initial -profile is quantified through the root-meansquare of the difference between the spatial profile of initial nucleation length ini c ( ) related to ini( ) and the perturbed one pert c ( ) related to ini( ) +̂( ) . (c-e) When the perturbation wavelength pert is smaller than the initial nucleation length ini c , it does not change its spatial profile, no matter the perturbation amplitude. Only the larger wavelength perturbations may influence the nucleation length.
10.1029/2021GL094901 9 of 11 To further demonstrate the potential of our theoretical framework, we consider a heterogeneous fault with a weakening rate profile ini( ) (see Figure 4a) that emerges from a multiscale distribution of asperities with a Hurst exponent = 0.7 (Ampuero et al., 2006). The nucleation length ini c ( ) can be computed from Equation 4. We superpose to the initial weakening rate profile ini , a unimodal perturbation ̂ of period pert , and amplitude pert , giving birth to a perturbed profile pert : When computing the nucleation length pert c ( ) associated to pert( ) , we observe that only perturbations whose wavelength is larger than the reference nucleation length ini c ( ) matter (see Figure 4b and 4e), and that the perturbation of critical length increases with both the wavelength pert and the amplitude pert . Furthermore, the small-scale perturbations ( pert ≲ hom c ) do not change the nucleation length pert c , whether its amplitude is small (Figure 4c) or large (Figure 4d). Note that if the initial nucleation length ini c locally drops to extremal values (see Figure 3d), very small-scale asperities may then influence the overall stability behavior of the fault.
Overall, our work provides then quantitative reasoning to assess which scale of asperities should be included in the modeling of complex faults and which can be averaged, when the frictional heterogeneities span over several length scales.

Conclusion
Nucleation processes along fault with differential frictional weakening are a collective phenomenon that may involve the progressive depinning of multiple asperities until a perturbation of size c is reached. Building on the theory of static friction and the physics of depinning, we proposed an analytical framework that allows to predict the effective critical length c for any spatial profile of . This framework has been successfully compared to dynamic simulations of Mode II friction and is directly tractable to nucleation of Mode I and Mode III fracture along weak interfaces. It provides clues to explain various nucleation scenarios observed in laboratory experiments and in nature as well as to derive scale-separation conditions assessing the influence of one asperity scale on the overall fault stability. Furthermore, it may provide ways to estimate the shear loading levels and the position at which nucleation occurs along more complex interfaces of known frictional properties, where all frictional quantities ( p, r , ) as well as the external loading ( n, ∞) fluctuate due to, for example, fault roughness (Cattania & Segall, 2021). The recent analogy drawn between earthquake nucleation for rate-and-state friction and rupture initiation along heterogeneous piecewise linear slip-weakening interfaces (Ray & Viesca, 2017;Viesca, 2016a) provides convincing ways to extend the proposed framework to rate-and-state friction laws. The generalization of our results to (a) nonlinear slip-dependent friction laws and (b) a three-dimensional setting is not straightforward, but one may adapt the approach proposed in this work to handle fault inhomogeneities to the energetic nucleation framework of Rice and Uenishi (2010). Yet, further experimental work is needed to assess the validity of our framework in predicting the influence of heterogeneities on the nucleation process. Dynamic rupture experiments performed on model microarchitectured faults in the laboratory (e.g., Berman et al., 2020) may constitute a critical test for our theory.

Data Availability Statement
Data regarding the figures of the main text are available on Zenodo (Lebihain et al., 2021).