Earthquake Nucleation Along Faults With Heterogeneous Weakening Rate.

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 L c 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.

Introduction

Understanding how interfaces fail is of utmost importance in fields ranging from earthquake physics to engineering fracture mechanics. For unstable frictional interfaces such as seismogenic faults, field observations (Bouchon et al., 2013; Kato et al., 2012) as well as laboratory experiments (Ben‐David & Fineberg, 2011; Dieterich, 1978; Latour et al., 2013; McLaskey, 2019; Ohnaka & Kuwahara, 1990) suggest one possible scenario where the onset of fault motion is characterized by the transition from quasistatic slip growth to dynamic rupture propagation (Passelègue et al., 2016; Svetlizky et al., 2016). The transition happens when a region of critical size of the fault is slipping. The knowledge of this nucleation length proves crucial since it allows to predict both the loading levels and the position at which the earthquake motion starts (Albertini et al., 2020; Ampuero et al., 2006; Uenishi & Rice, 2003).

Previous theoretical works linked to the frictional properties of the fault for linear slip‐dependent (Campillo & Ionescu, 1997; Dascalu et al., 2000; Uenishi & Rice, 2003) and more complex rate‐and‐state (Aldam et al., 2017; Brener et al., 2018; Ruina, 1983; Rubin & Ampuero, 2005; Viesca, 2016b) friction laws along homogeneous faults. Yet, frictional properties are expected to significantly vary along the fault plane and 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 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.

Materials and Methods

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 , 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 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 to its residual one with a weakening rate that describes the material brittleness/ductility. Variations of the frictional properties 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.

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 of the interface goes from its peak value to the residual one when the local slip reaches its critical value , defining a weakening rate . , , 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 . (c) Slip growth develops quasistatically without any significant velocity burst, until one slip patch reaches a critical length that leads to the dynamic rupture of the whole interface (see Movie S1).

As the macroscopic loading grows, it locally exceeds the friction , 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, the shear wave velocity, ( and being the shear modulus and the Poisson's ratio, respectively), and is a linear operator. In Equation 1, the term , 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 is low. As the loading is further increased, they grow quasistatically and coalesce into larger slipping regions. This initial nucleation stage of duration proves rather quiescent since no major velocity burst is observed. Yet, at , 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 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 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 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 and the residual friction are uniform (Figure 2a). To make any parallel to Mode I fracture easier and without any loss of generality, we set (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 between two extremal values . We set the nucleation length of the reference homogeneous material with uniform as the adimensionalizing length of the system (Uenishi & Rice, 2003). In the following, we consider , , , and . 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.

Figure 2

Measuring the nucleation length of a heterogeneous fault: (a) and 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 between two extremal values . (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) 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 (in linear‐log space), which provides an accurate measurement of from the growth rate . See Movie S2, Movie S3, and Figure S2 in the Supporting Information S1 for comparison with the homogeneous case of Uenishi and Rice (2003).

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 (), where the stress locally exceeds the frictional resistance (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 where a portion of the interface is slipping, while (b) the second involves an unstable dynamic crack propagation for 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 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 ( in Figure 2).

To quantify the influence of spatial variations of on , 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 , as previously observed in laboratory experiments of earthquake nucleation between two polycarbonates blocks (Latour et al., 2013). The nucleation length may then be estimated from the maximal growth rate , as the length where the patch expansion is at its strongest. The validity of our heuristic approach is assessed on homogeneous interfaces for which is known a priori (Favreau et al., 1999; Uenishi & Rice, 2003). We use it to numerically estimate the critical length of heterogeneous interfaces with a precision, corresponding to the error observed for the homogeneous interface of known (see Section 2.1 in Supporting Information S1).

Results and Discussion

Theoretical Model for Nucleation Length Predictions

Here, we propose a way to determine the critical length analytically building on both the theory of static friction and the physics of depinning. For a velocity perturbation centered in with a support of size , Equation 2 becomes: where is the linear operator introduced in Dascalu et al. (2000) and Uenishi and Rice (2003), and 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 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 up to the second order and solve the eigenproblem , where (see Section 3 in Supporting Information S1). The nucleation length is the solution of the transcendental Equation 4, which encompasses the main novelty of the paper.

In Equation 4, denotes the eigenmode associated to the eigenvalue of the homogeneous eigenproblem (Dascalu et al., 2000; Uenishi & Rice, 2003). The first two terms of Equation 4 represent the heterogeneity contributions up to the first order. The value of the critical length at a position involves spatial variations of on a scale potentially larger than the heterogeneity size . This collective yet heterogeneous behavior in the earthquake nucleation cannot be grasped by the homogeneous theory (Albertini et al., 2020; Uenishi & Rice, 2003). The third term corresponds to the second‐order contributions up to a critical mode . 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., ). 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 of Equation 4 to the numerically estimated critical lengths 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).

Figure 3

(a) The length of the critical instability nucleus measured from the dynamic simulations is compared to the theoretical prediction of Equation 4 at the fault center for a broad range of characteristic scale of the asperities (320 simulations). Red region: error on . (b) may strongly differ from the local nucleation length of Uenishi and Rice (2003) at . 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 set by the average frictional properties, the effective nucleation length of Equation 4 follows Uenishi and Rice (2003)'s predictions ; (d) when is smaller than yet larger than the minimal nucleation length set by the most brittle defect, departs significantly from but can be locally controlled by the extrema of the weakening rate distribution . Inset: instability birth at ; (e) when is smaller than both and , the nucleation behavior is homogenized and the nucleation length is comparable to that set by the average frictional properties (see Movies [Link], [Link], [Link]).

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 along the interface (see black solid lines in Figures 3c–3e). In our simulations, the effective nucleation length corresponds to the one predicted at due to the peaked nature of the loading. In more realistic cases, the position of the earthquake nucleus (and the associated nucleation length ) will depend not only (a) on heterogeneities of peak strength , 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.

Figure 4

Influence of a modal perturbation on the stability of a heterogeneous fault: (a) an initial weakening rate profile consisting of the superposition of multiple spatial modes is considered. The spatial variations of give birth to an effective profile of nucleation length following Equation 4. (b) The influence of the superposition of a unimodal perturbation of period and amplitude to the initial ‐profile is quantified through the root‐mean‐square of the difference between the spatial profile of initial nucleation length related to and the perturbed one related to . (c–e) When the perturbation wavelength is smaller than the initial nucleation length , it does not change its spatial profile, no matter the perturbation amplitude. Only the larger wavelength perturbations may influence the nucleation length.

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 , and the scale set by the weakest defect along the fault .

Local regime: when and , 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 , which can be distributed above () or below ()

Extremal regime: when and , a critical nucleation patch of size may develop within a single brittle asperity of size , where the weakening rate reaches its maximal value . 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 ). Along natural faults, these small ruptures may be arrested by local barriers of strength , 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 far away from them, and that, in contrast, and (b) other weak spots may not be critical if not brittle enough that is., (see the spatial evolution of in black line in Figure 3d and Section 4 in Supporting Information S1)

Homogenized regime: when and , no critical slip patch can develop within a single asperity. Nucleation occurs after the collective depinning of multiple asperities (Dublanchet, 2018; Perfettini et al., 2003; Ray & Viesca, 2019) with dynamics similar to that of homogeneous nucleation. The critical length fluctuates around its homogenized value set by the averaged frictional properties (see Figure 3e) and can then be studied within the homogeneous nucleation theory of Uenishi and Rice (2003). Note that Equation 4 fully captures the fluctuations of around , which may grow as the intensity of weakening rate fluctuations increases

We argue here that all three instability regimes could occur along natural faults depending on their size, geometry, maturity, and lithology. But the homogenized regime proves to be of major importance for geophysical applications. Indeed, heterogeneous fracture is often described as a critical phenomenon controlled by the weakest defect, thus ruling out its study within a homogeneous framework. Yet, our results suggest that under the scale‐separation condition , the stability behavior of a heterogeneous fault can be studied within the homogeneous framework of Favreau et al. (1999) and Uenishi and Rice (2003) with . Moreover, the existence of the homogenized regime may account for the relative reproductivity of laboratory experiments where sample roughness is often imposed and kept relatively smooth and further justifies their relevance in the modeling of natural faults.

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.

To further demonstrate the potential of our theoretical framework, we consider a heterogeneous fault with a weakening rate profile (see Figure 4a) that emerges from a multiscale distribution of asperities with a Hurst exponent (Ampuero et al., 2006). The nucleation length can be computed from Equation 4. We superpose to the initial weakening rate profile , a unimodal perturbation of period , and amplitude , giving birth to a perturbed profile :

When computing the nucleation length associated to , we observe that only perturbations whose wavelength is larger than the reference nucleation length matter (see Figure 4b and 4e), and that the perturbation of critical length increases with both the wavelength and the amplitude . Furthermore, the small‐scale perturbations () do not change the nucleation length , whether its amplitude is small (Figure 4c) or large (Figure 4d). Note that if the initial nucleation length 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 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 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 as well as the external loading 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.

Supporting Information S1

Click here for additional data file.

Movie S1

Click here for additional data file.

Movie S2

Click here for additional data file.

Movie S3

Click here for additional data file.

Movie S4

Click here for additional data file.

Movie S5

Click here for additional data file.

Movie S6

Click here for additional data file.