Stability of pulse-like earthquake ruptures

TL;DR

This paper investigates the stability of pulse-like earthquake ruptures driven by thermal pressurisation, revealing their inherent instability and implications for fault strength and rupture complexity through simulations and theoretical modeling.

Contribution

It introduces a dynamic elastodynamic model showing that steady-state slip pulses driven by thermal pressurisation are unstable, providing insights into rupture evolution and fault strength.

Findings

Steady-state pulse-like ruptures are unstable under thermal pressurisation.

Positive stress perturbations lead to rupture growth and transition to other modes.

Unstable pulses influence minimum stress conditions for sustainable ruptures.

Abstract

Pulse-like ruptures arise spontaneously in many elastodynamic rupture simulations and seem to be the dominant rupture mode along crustal faults. Pulse-like ruptures propagating under steady-state conditions can be efficiently analysed theoretically, but it remains unclear how they can arise and how they evolve if perturbed. Using thermal pressurisation as a representative constitutive law, we conduct elastodynamic simulations of pulse-like ruptures and determine the spatio-temporal evolution of slip, slip rate and pulse width perturbations induced by infinitesimal perturbations in background stress. These simulations indicate that steady-state pulses driven by thermal pressurisation are unstable. If the initial stress perturbation is negative, ruptures stop; conversely, if the perturbation is positive, ruptures grow and transition to either self-similar pulses (at low background stress) or expanding cracks (at elevated background stress). Based on a dynamic dislocation model, we develop an elastodynamic equation of motion for slip pulses, and demonstrate that steady-state slip pulses are unstable if their accrued slip $b$ is a decreasing function of the uniform background stress $\tau_\mathrm{b}$. This condition is satisfied by slip pulses driven by thermal pressurisation. The equation of motion also predicts quantitatively the growth rate of perturbations, and provides a generic tool to analyse the propagation of slip pulses. The unstable character of steady-state slip pulses implies that this rupture mode is a key one determining the minimum stress conditions for sustainable ruptures along faults, i.e., their ``strength''. Furthermore, slip pulse instabilities can produce a remarkable complexity of rupture dynamics, even under uniform background stress conditions and material properties.

Figures

Equations

$$\tau (X)={\tau }_{\mathrm{b}}-\frac{\bar{\mu }}{2\pi v_{\mathrm{r}}}{\int }_0^L\frac{V(\xi )}{X-\xi }d\xi ,$$

$$V(X)=v_{\mathrm{r}}\frac{d\delta }{dX},$$

$$k_{\mathrm{L}}=-\frac{4}{\pi \sqrt{L}}{\int }_0^L\sqrt{\frac{X}{L-X}}\frac{d\tau }{dX}dX.$$

$${\tau }_{\mathrm{f}}=f{\sigma }^{\prime }=f\times ({\sigma }_{\mathrm{n}}-p),$$

$$\frac{\partial p}{\partial t}$$

$$\frac{\partial \Theta }{\partial t}$$

$$\dot{\gamma }(y,t)=\frac{V(t)}{h}{e}^{-\pi y^2/h^2},$$

$${\tau }_{\mathrm{f}}={\tau }_0-\frac{1}{{\delta }_{\mathrm{c}}}{\int }_0^t{\tau }_{\mathrm{f}}(t^{\prime })V(t^{\prime })\mathcal{K}\left(\frac{t-t^{\prime }}{T^{\ast }};\frac{{\alpha }_{\mathrm{hy}}}{{\alpha }_{\mathrm{th}}}\right)d{t}^{\prime },$$

$${\delta }_{\mathrm{c}}=\frac{\rho c}{f\Lambda }h\text{and}{T}^{\ast }=\frac{h^2}{4\alpha },$$

$$V^{\ast }={\delta }_{\mathrm{c}}/T^{\ast }\text{and}{L}^{\ast }=\mu {\delta }_{\mathrm{c}}/{\tau }_0.$$

$$h_{\mathrm{dyna}}=\frac{\mu }{{\tau }_0}\frac{\rho c}{f\Lambda }\frac{4\alpha }{c_{\mathrm{s}}},$$

$$\tau (x,t)={\tau }_{\mathrm{b}}(x)-\frac{\mu }{2c_{\mathrm{s}}}V(x,t)+\varphi [V],$$

$$\Delta \tau (x,t)=\Delta {\tau }_{\mathrm{b}}(x)-\frac{\mu }{2c_{\mathrm{s}}}\Delta V(x,t)+\varphi [\Delta V],$$

$$\Delta {\tau }_{\mathrm{f}}(x,t)$$

$$+\Delta {\tau }_{\mathrm{f}}(x,t^{\prime })\Delta V(x,t^{\prime }))\mathcal{K}\left(\frac{t-t^{\prime }}{T^{\ast }};\frac{{\alpha }_{\mathrm{hy}}}{{\alpha }_{\mathrm{th}}}\right)d{t}^{\prime },$$

$$\tau (x,t)={\tau }_{\mathrm{b}}-\frac{\mu }{2c_{\mathrm{s}}}V(x,t)+\varphi (x,t),$$

$$\varphi (x,t)=\frac{\mu }{2\pi }\frac{\partial }{\partial x}{\int }_{-\infty }^t{\int }_{-\infty }^{+\infty }M\left(\frac{x-x^{\prime }}{c_{\mathrm{s}}(t-t^{\prime })}\right)\frac{\partial \delta /\partial x^{\prime }(x^{\prime },t^{\prime })}{t-t^{\prime }}d{x}^{\prime }d{t}^{\prime }.$$

$$M(u)=\mathcal{H}(1-u^2)\sqrt{1-u^2},$$

$$|x-\xi (t)|\gg L(t):$$

$$\frac{\partial \delta }{\partial x}(x,t)=-b(\theta (x)){\delta }_{\mathrm{Dirac}}(\xi (t)-x)$$

$$+\frac{db(\theta (x))}{dx}\mathcal{H}(\xi (t)-x)$$

$$\tau (x,t)-{\tau }_{\mathrm{b}}=\varphi (x,t)-\frac{\mu }{2c_{\mathrm{s}}}{v}_{\mathrm{r}}(t)b(\theta (x)){\delta }_{\mathrm{Dirac}}(\xi (t)-x)$$

$$\varphi (x,t)$$

$${\varphi }_{\mathrm{Dirac}}(x,t)$$

$${\varphi }_H(x,t)$$

$$\bar{u}\equiv \frac{x-\xi (t^{\prime })}{c_{\mathrm{s}}(t-t^{\prime })}.$$

$$|x-\xi (t)|\gg L(t):\tau (x,t)-{\tau }_{\mathrm{b}}=\varphi (x,t).$$

$$X=\xi (t)-x.$$

$$\tau (x,t)-{\tau }_{\mathrm{b}}={\varphi }_{\mathrm{ss}}(X;v_{\mathrm{r}},b)\equiv -\frac{\mu \sqrt{1-v_{\mathrm{r}}^2/c_{\mathrm{s}}^2}}{2\pi }\frac{b}{X}.$$

$$L\ll L_{\mathrm{out}}.$$

$${\varphi }_{\mathrm{Dirac}}(X,t)$$

$$+\frac{\mu }{4\pi c_{\mathrm{s}}b(t)}\frac{d}{dt}\left[b^2(t)\frac{v_{\mathrm{r}}(t)/c_{\mathrm{s}}}{\sqrt{1-v_{\mathrm{r}}^2(t)/c_{\mathrm{s}}^2}}\right]\ln \left|\frac{X}{L_{\mathrm{out}}}\right|.$$

$${\varphi }_H(x,t)={\int }_{-\infty }^t(\partial {\varphi }_H(x,t)/\partial t)dt,$$

$$\frac{\partial {\varphi }_H(x,t)}{\partial t}=\frac{\mu }{2\pi c_{\mathrm{s}}}{\int }_{-\infty }^t\frac{dM\left(\bar{u}\right)}{d\bar{u}}\frac{db/d{t}^{\prime }}{(t-t^{\prime }{)}^2}d{t}^{\prime },$$

Full text

\setkeys

Gindraft=false \authorrunningheadBRANTUT ET AL. \titlerunningheadSLIP PULSE STABIILITY \authoraddrN. Brantut, Department of Earth Sciences, University College London, Gower Street, London WC1E 6BT, UK. (n.brantut@ucl.ac.uk)

11affiliationtext: Department of Earth Sciences, University College London, London, UK.22affiliationtext: Department of Civil and Resource Engineering, Dalhousie University, Halifax, Canada.33affiliationtext: Disaster Prevention Research Institute, Kyoto University, Uji, Japan.

{article}

1 Introduction

The propagation of earthquakes is generally classified into two main modes: crack-like ruptures, where fault slip occurs throughout the duration of propagation, and pulse-like ruptures, where only a small portion of the fault inside a ruptured area is sliding at a given time during rupture. The observation that local slip duration is often much shorter than the time required for stopping phases to propagate from fault boundaries led Heaton (1990) to suggest that most crustal earthquakes may propagate as pulse-like ruptures. A number of detailed kinematic and dynamic inversions of earthquake slip (Wald and Heaton, 1994; Beroza and Ellsworth, 1996; Olsen et al., 1997; Day et al., 1998; Galetzka et al., 2015) have confirmed the pulse-like nature of large crustal earthquakes, highlighting the importance of this rupture mode in the physics of faults.

The physical origin and dynamics of pulse-like ruptures have been studied extensively in theoretical models. Slip events have been shown to propagate as narrow, self-similar slip pulses in simplified discrete spring block models (e.g., Carlson and Langer, 1989; Elbanna and Heaton, 2012). Fully dynamic rupture simulations have revealed the key role of velocity-weakening friction (e.g. Heaton, 1990; Cochard and Madariaga, 1994; Perrin et al., 1995; Zheng and Rice, 1998) and boundary conditions (e.g. Johnson, 1990, 1992) in the spontaneous generation of slip pulses. Specifically, elastodynamic simulations with velocity-dependent friction show that the existence and evolution of the dynamic pulse-like ruptures are strongly controlled by both the ambient background stress and the nucleation conditions on the fault (Zheng and Rice, 1998; Gabriel et al., 2012). For a given nucleation condition, an increase in background stress results in a sequential transition from arresting pulses to growing pulses, and then growing crack-like ruptures. Therefore, the mode of rupture and its evolution are the signature of the background stress acting on the fault prior to the earthquake. Of critical importance here is the stress level at the transition from arresting to growing pulses, which provides the threshold below which sustained fault slip is precluded. The rupture mode at the transition is that of a “steady-state” slip pulse, for which the tip and tail of the slipping patch propagate at the same speed. These steady-state solutions are therefore key to understand the stress level required for earthquake propagation and the dynamics of faults.

Steady-state solutions of the elastodynamic fault problem can be obtained using analytical or simple numerical methods, so that they can be studied efficiently without resorting to computationally expensive numerical treatment. Several pulse-like rupture solutions have been obtained for simple models of faults with a constant or slip-dependent friction law (e.g. Broberg, 1978; Freund, 1979; Rice et al., 2005; Dunham and Archuleta, 2005), but without specific regards to the processes allowing for strength recovery and “healing” (i.e., cessation of slip) at the tail of the pulse. Steady-state pulse solutions fully consistent with both elastodynamics and a specific friction law have been determined by Perrin et al. (1995) in the context of rate-and-state friction, and more recently by Garagash (2012) and Platt et al. (2015) in the context of dynamic weakening by thermal (or chemical) pressurisation of pore fluids within the fault zone. These solutions provide unique insights into the relationships between rupture properties, such as pulse width or rupture velocity at a given background stress, and key parameters of the friction law, such as rate-and-state parameters (Perrin et al., 1995) or thermo-hydraulic properties of the fault core (Garagash, 2012; Platt et al., 2015). Despite the (relative) simplicity and efficiency of those steady state pulse solutions, it remains to be confirmed how they can be generated and how they evolve in response to perturbations in loading conditions or frictional properties. In other words, the key question here is to determine how self-consistent steady-state solutions (i.e., satisfying elastodynamics and all the features of a specific friction law) can be compared to possibly transient rupture dynamics observed on natural faults.

Regarding this issue, the numerical simulations provided by Gabriel et al. (2012) and Brener et al. (2018) using velocity-dependent friction, or by Noda et al. (2009) in the context of weakening by thermal pressurisation of pore fluids within the fault, seem to indicate that such steady-state solutions are not stable: they either grow (to form self-similar pulses) or decay and stop. The goal of this paper is to analyse in detail how steady-state pulses respond to perturbations and to determine a clear stability condition depending on the characteristics of the friction law. Building on the work by Garagash (2012), we examine specifically the case of pulses driven by thermal pressurisation of pore fluids, and first solve the nonlinear perturbation problem numerically (Section 2). We then examine more generally the conditions under which stable pulses can exist based on an approximate equation of motion for moving dislocations (Section 3). The significance of steady-state pulse solutions and some implications for the dynamics of earthquakes are examined in Section 4.

2 Slip pulses driven by thermal pressurisation of pore fluids

In this Section, we present a detailed analysis of the evolution of pulses driven by thermal pressurisation. We choose to focus specifically on thermal pressurisation as the governing process by which faults weaken (and restrengthen), since it has a firm physical background, and has been shown to be consistent with a number of seismological observations (Rice, 2006; Viesca and Garagash, 2015). Beyond this specific choice for the fault constitutive behaviour, we stress that the method of analysis developed here is quite general and can be used to include other friction laws.

We first briefly summarise the results of Garagash (2012) regarding steady-state solutions, and perform a stability analysis by solving for the evolution of perturbations from the steady-state solution.

2.1 Model and steady-state solution

2.1.1 Elastodynamics of steadily-propagating pulse

We consider a planar fault embedded in an infinite, homogeneous elastic medium of shear modulus $\mu$. The fault is assumed to be of infinite extent in one of its planar dimensions, so that we restrict our attention to a two-dimensional problem. The fault is loaded by a uniform background shear stress, denoted $\tau_{\mathrm{b}}$. For simplicity, we assume that the loading is in mode III (out of plane) geometry. Fault slip is assumed to occur over a patch of finite length $L$, which propagates at a constant speed $v_{\mathrm{r}}$ along spatial coordinate $x$, as shown in Figure 1. Under steady-state conditions (i.e., constant rupture speed), it is convenient to introduce a reference frame $(X,y)$ that moves with the rupture tip, so that the shear stress $\tau$ and slip rate $V$ along the fault (in the plane $y=0$) are functions of the coordinate $X=v_{\mathrm{r}}t-x$ only. The elastodynamic equilibrium requires that (e.g. Weertman, 1969)

[TABLE]

where $\bar{\mu}$ is an apparent shear modulus given by $\bar{\mu}=\mu\times F(v_{\mathrm{r}}/c_{\mathrm{s}})$. The function $F$ of the ratio of rupture speed and shear wave speed $c_{\mathrm{s}}$ is equal to $F(v_{\mathrm{r}}/c_{\mathrm{s}})=\sqrt{1-v_{\mathrm{r}}^{2}/c_{\mathrm{s}}^{2}}$ ((e.g. Rice, 1980)), so that the apparent modulus approaches zero as the rupture speed approaches the shear wave speed. In Equation (1), it is understood that the slip rate is given by

[TABLE]

where $\delta$ is the slip.

In the slipping part of the fault ($0\leq X\leq L$), the stress $\tau(X)$ must be equal to the fault strength $\tau_{\mathrm{f}}$, which is given by a constitutive law (see below). Furthermore, at the tail of the pulse ($X\geq L$), we need to ensure that the strength remains higher than the elastic stress $\tau(X)$ (otherwise slip would continue, which would be in contradiction with the pulse width being equal to $L$). Garagash (2012) determined that the stress gradient at the tail of the pulse is singular, of the form $d\tau/dX\propto k_{\mathrm{L}}/\sqrt{X-L}$, where

[TABLE]

The gradient in fault strength remains continuous, so that the condition for cessation of slip $\tau(X)\leq\tau_{\mathrm{f}}$ imposes that the elastic stress gradient remains bounded, i.e., $k_{\mathrm{L}}=0$. This equality ensures the consistency of the assumption that slip only occurs where $\tau_{\mathrm{f}}=\tau(X)$, and provides a constraint on the pulse length $L$ (which would otherwise be a free parameter of the problem).

2.1.2 Fault strength

The strength of the fault $\tau_{\mathrm{f}}$ is assumed to be governed by a friction law:

[TABLE]

where $f$ is a friction coefficient and $\sigma^{\prime}$ is the Terzaghi effective stress, equal to the difference between the fault normal stress $\sigma_{\mathrm{n}}$ and the pore fluid pressure $p$ inside the fault core. Here, we assume a constant friction coefficient throughout the slip process. Since we are primarily concerned here with dynamic slip, the constant value of $f$ should be representative of the high velocity “dry” friction coefficient, which is typically of the order of $0.1$ (e.g. Di Toro et al., 2011). The pore fluid pressure $p$ is governed by the competition between fluid diffusion and thermal expansion due to shear heating. The fluid pressure evolution is coupled to the temperature $\Theta$ through the following Equations (e.g. Rice, 2006):

[TABLE]

where $\Lambda$ is a thermo-poro-elastic coupling factor expressing the increase in fluid pressure per unit increase in temperature, $\alpha_{\mathrm{hy}}$ and $\alpha_{\mathrm{th}}$ are hydraulic and thermal diffusivities, respectively, $\dot{\gamma}$ is the shear strain rate, and $\rho c$ is the heat capacity of the fault rock. In a fault core of finite width, pore pressure is not homogeneous across the fault, which can lead to shear strain localisation (Rice et al., 2014; Platt et al., 2014). Here, we do not explicitly account for this effect, which requires the introduction of further parameters such as rate-hardening properties of the sheared gouge. We follow Garagash (2012) and consider a Gaussian shear strain rate distribution across the fault with a characteristic width $h$,

[TABLE]

and use the pore pressure at the center of the fault (where it is maximum) to compute the strength in Equation (4).

Equations (5) and (6) can be solved to arrive at the integral representation for fault strength (Rice, 2006) given here in the form of (Garagash, 2012):

[TABLE]

where $\tau_{0}=f\sigma^{\prime}_{0}$ is the initial strength of the fault (at $p(y=0,t=0)=p_{0}$), and

[TABLE]

where $\alpha=(\sqrt{\alpha_{\mathrm{th}}}+\sqrt{\alpha_{\mathrm{hy}}})^{2}$, are characteristic slip weakening distance and diffusion time, respectively. The convolution kernel $\mathcal{K}$ is given in Appendix A.

2.1.3 Steady-state pulse solution

For a given background stress $\tau_{\mathrm{b}}/\tau_{0}$ and diffusivity ratio $\alpha_{\mathrm{hy}}/\alpha_{\mathrm{th}}$, Equations (1) and (8), under the condition (3), have been solved for slip rate $V$ and strength $\tau_{\mathrm{f}}$ by Garagash (2012). Here, we reproduce these computations using a more efficient quadrature method given by Viesca and Garagash (2018) (see Appendix B.1 for more details). To keep the solutions as general as possible, we normalise the stresses by $\tau_{0}$, slip by $\delta_{\mathrm{c}}$, time by $T^{*}$, and slip rate and distance by

[TABLE]

In the determination of the solution, we constrain not only the distribution of stress and slip rate along the pulse, but also its duration $T/T^{*}$ and length $L/L^{*}$. The rupture speed is given by $v_{\mathrm{r}}=L/T$, so that its ratio relative to the shear wave speed is $v_{\mathrm{r}}/c_{\mathrm{s}}=(L/L^{*})(T^{*}/T)/(c_{\mathrm{s}}T^{*}/L^{*})$. Following Garagash (2012, Section 7.1), we define a characteristic thickness $h_{\mathrm{dyna}}$ such a that $h/h_{\mathrm{dyna}}=c_{\mathrm{s}}T^{*}/L^{*}$, i.e.,

[TABLE]

so that constraining the fault core thickness through the ratio $h/h_{\mathrm{dyna}}$ implies that the rupture velocity $v_{\mathrm{r}}/c_{\mathrm{s}}$ is also constrained ($v_{\mathrm{r}}/c_{\mathrm{s}}=(L/L^{*})(T^{*}/T)(h_{\mathrm{dyna}}/h)$).

A representative example is shown in Figure 2, where we chose $\tau_{\mathrm{b}}/\tau_{0}=0.7$ and $\alpha_{\mathrm{hy}}/\alpha_{\mathrm{th}}=1$. For completeness, we also show the evolution of stress $\tau$ and strength $\tau_{\mathrm{f}}$ outside the pulse, and we indeed observe that $\tau<\tau_{\mathrm{f}}$ behind the tail. This pulse is therefore fully consistent with elastodynamics and the fault constitutive law.

Some key properties of steady-state pulses driven by thermal pressurisation can be determined from a systematic exploration of the numerical solutions. Of particular interest here are the pulse width ($L/L^{*}$), total slip ($b/\delta_{\mathrm{c}}$) and rupture speed ($v_{\mathrm{r}}/c_{\mathrm{s}}$), which are shown in Figure 3 as a function of the background stress ($\tau_{\mathrm{b}}/\tau_{0}$). In all these plots, we chose again $\alpha_{\mathrm{hy}}/\alpha_{\mathrm{th}}=1$, knowing that this parameter has only a minor quantitative effect on the results (Garagash, 2012).

With increasing background stress, the slip and rupture speed decrease, while the pulse width increases. However, the relationships depicted in Figure 3 have been derived from independent steady-state solutions, and therefore may not correspond to the actual evolution of the width, slip and rupture speed of a single pulse propagating along a fault with varying background stresses, pore pressure, or fault constitutive parameters (friction, etc). In order to compute such an evolution, and determine whether a given steady-state solution is stable against perturbations in background stress, we need to compute the full elastodynamic solution for a propagating pulse in a perturbed stress state.

2.2 Elastodynamic stability analysis: Method

The elastodynamic stress equilibrium can be expressed as

[TABLE]

where $\phi[V]$ is a linear functional of slip rate that corresponds to the stress redistribution due to slip and elastic waves. In Equation (12), an explicit space dependency has been written for the background stress, $\tau_{\mathrm{b}}(x)$, to account for the introduction of local perturbations. Direct solutions of Equation (12) can be obtained numerically, but require somewhat arbitrary rupture initiation conditions, which would be incompatible with our objective of studying small perturbations around a steadily propagating rupture, regardless of how this rupture originated. Here, we circumvent the rupture nucleation problem and only solve for stress and slip rate perturbations from a preexisting steady-state pulse solution.

Let us denote $\tau_{\mathrm{ss}}(x,t)$, $\delta_{\mathrm{ss}}(x,t)$ and $V_{\mathrm{ss}}(x,t)$ the stress, slip and slip rate associated with a steady-state pulse propagating along a fault under a uniform background stress $\tau_{\mathrm{b,ss}}$. By construction, $\tau_{\mathrm{ss}}$, $\delta_{\mathrm{ss}}$ and $V_{\mathrm{ss}}$ are solutions of Equation (12) with $\tau_{\mathrm{b}}=\tau_{\mathrm{b,ss}}$. Now if we introduce a perturbation in background stress $\Delta\tau_{\mathrm{b}}(x)$, the resulting perturbations $\Delta\tau(x,t)$, $\Delta\delta(x,t)$ and $\Delta V(x,t)$ in stress, slip and slip rate, respectively, satisfy

[TABLE]

where we made use of the linearity of the functional $\phi$. The strength evolution due to thermal pressurisation is given in Equation (8), which is rewritten in terms of strength and slip rate perturbations as

[TABLE]

where $\tau_{\mathrm{f,ss}}(x,t)$ is the strength along the steady-state pulse and $\Delta\tau_{\mathrm{f}}$ is the strength perturbation. The governing Equation (14) for the strength perturbation is not linear, and therefore requires the specific knowledge of the steady-state strength and slip rate profiles, $\tau_{\mathrm{f,ss}}(x,t)$ and $V_{\mathrm{ss}}(x,t)$. These profiles correspond to the solutions of the steady-state problem stated in the previous Section.

Our solution strategy therefore consists in first solving a steady-state problem (see previous Section), and then solving the full elastodynamic problem for perturbations to this solution arising from variations in background stress. In practice, we use the spectral boundary integral method of Perrin et al. (1995); Lapusta et al. (2000); Noda and Lapusta (2010) to compute the dynamic stress distribution functional $\phi[\Delta V]$ (see Appendix B.2), and a predictor-corrector method for time integration. The details of the algorithm are given in Appendix B.3.

2.3 Elastodynamic stability analysis: Results

Two representative examples of slip pulse propagating across either a positive or a negative perturbation in background stress are shown in Figure 4. The initial background stress is $\tau_{\mathrm{b}}/\tau_{0}=0.7$ and the diffusivity ratio is $1$. In both cases, the perturbation was a half-sine of $(1/3)L/L^{*}$ in width and $\Delta\tau_{\mathrm{b}}/\tau_{0}=10^{-2}$ in amplitude. When crossing a negative perturbation (Figure 4a,b), the slip pulse continues to propagate over a distance of the order of $10L^{*}$ while both the dynamic stress drop and slip rate progressively reduce, until rupture arrests. Conversely, a positive perturbation (Figure 4c,d) amplifies the dynamic stress drop and slip rate, and also results in a progressive increase in pulse width.

The evolution of the pulse shape is best observed in the coordinate system that moves with the pulse tip at its reference speed $v_{\mathrm{r,ss}}$. Figure 5 shows contours of slip rate in the transformed coordinate system \big{(}(x-v_{\mathrm{r,ss}}t)/L^{*},t/T^{*}\big{)} for the two simulations presented in Figure 4. When the perturbation is negative (Figure 5a), the pulse width reduction is initially driven by an acceleration of the trailing edge (healing front), and subsequently by a deceleration of the tip. The acceleration of the healing front initiates when the shear wave emitted from the pulse tip at the location of the perturbation reaches the trailing edge of the pulse. The overall pulse width reduces in a nonlinear manner over time, and the pulse arrests abruptly. When the perturbation is positive (Figure 5b), an acceleration of the pulse tip is first observed, followed by a deceleration of the trailing edge. The pulse tip speed gradually approaches the shear wave speed. After a critical time of the order of $\sim 20T^{*}$, the trailing edge accelerates again and further propagates at a speed greater than the

initial $v_{\mathrm{r,ss}}$, but less than the tip speed. The slip pulse then becomes expanding.

The time evolution of perturbations in normalised pulse width, peak slip rate and slip is given in Figure 6. Regardless of the sign of the perturbation in background stress, all the perturbed quantities appear to grow exponentially with time since the perturbation onset $t_{\mathrm{pert}}$ at which the pulse tip enters the perturbed region, until either the complete arrest of the pulse (Figure 6a) or the transition to an expanding pulse (Figure 6b, at $(t-t_{\mathrm{pert}})/T^{*}\gtrsim 25$). The amplitude of all the normalised variables is initially of the same order of magnitude as the perturbation in background stress, and all grow at approximately the same exponential rate. As shown in Figure 14, the amplitude of the perturbation impacts only the initial jump in the perturbed variables but does not modify the growth rate itself.

The growth rate of slip perturbations following a negative stress perturbation is explored as a function of reference background stress in Figure 7. Reasonably accurate simulations can only be performed for $\tau_{\mathrm{b}}/\tau_{0}\gtrsim 0.4$, because at lower stress levels the reference rupture speeds becomes too close to the shear wave speed (see Figure 3c). A clear trend of increasing growth rate with increasing reference shear stress is observed. This trend is not linear (Figure 7b): the growth rate approaches zero at low stress, and increases dramatically at high stress. In reference to Figure 3c, the overall trend implies that slow pulses with little net slip arrest more rapidly than faster ones.

At elevated background stress ($\tau_{\mathrm{b}}/\tau_{0}\gtrsim 0.79$), the pulse response to positive perturbations is qualitatively different from that shown in Figures 4c,d and 5b. Figure 8 shows a series of snapshots of shear stress and slip rate for a pulse propagating under a background stress $\tau_{\mathrm{b}}/\tau_{0}=0.9$ and perturbed at $x=0$ with a half-sine of width $(1/3)L/L^{*}$ and amplitude $\Delta\tau_{\mathrm{b}}/\tau_{0}=10^{-3}$. The tip of the pulse accelerates, and the tail decelerates and starts propagating in the negative $x$ direction, leaving an expanding region of nonzero slip rates across the crack line. The pulse-like rupture effectively transitions to a crack-like rupture.

At intermediate background stress ($\tau_{\mathrm{b}}/\tau_{0}\sim 0.78$ for $\alpha_{\mathrm{hy}}/\alpha_{\mathrm{th}}=1$ and $h/h_{\mathrm{dyna}}=1$), a positive stress perturbation produces a complex rupture pattern, shown in Figure 9. The slip pulse initially transitions to a self-similar expanding pulse. At later times, a new, crack-like, rupture appears near the location where this transition occurred (Figure 9a). A plot of the shear stress and strength profiles (Figure 9b) reveals that the secondary nucleation is driven by the combination of a reduced strength in the wake of the pulse (although offset by strength recovery driven by pore fluid diffusion), and an increased backstress due to the expansion of the pulse. The net slip due to the expanding pulse increases approximately linearly with increasing propagation distance, so that the shear stress around the transition point from steady-state to expanding pulse is expected to increase logarithmically with time, and secondary nucleation ensues.

The process by which secondary nucleation might proceed is illustrated in detail in Figure 10, which shows (a) the stress profiles, and (b) the maximum stress perturbation behind the pulse (as well as the net slip perturbation) as a function of time for a simulation with positive stress perturbation and an initial background stress $\tau_{\mathrm{b}}/\tau_{0}=0.77$. Although secondary nucleation was not observed within the time frame of that simulation, a clear progressive increase in stress is observed around the location of the transition from steady-state to expanding pulse (see stress profiles inside the box in Figure 10(a)). A similar mechanism for the secondary rupture nucleation in the wake of expanding primary pulse, pulse or crack depending on the background stress level, have been described by Gabriel et al. (2012) for a fault with a velocity-weakening friction. This increase slows with increasing time and propagation distance (Figure 10(b)), but does not stabilise. This logarithmic increase in stress is expected if the net slip behind the pulse increases linearly with propagation distance, which seems to be the case here.

In summary, the numerical results presented above indicate that pulse-like ruptures driven by thermal pressurisation of pore fluids are unstable to infinitesimal perturbations. The growth of slip rate, slip and pulse width perturbations is initially exponential, and of the same sign as the initial stress perturbation. When that perturbation is negative, the slip pulse eventually stops, and does so abruptly. When the perturbation is positive, depending on the initial uniform background stress, the slip pulse grows and transitions to a self similar expanding pulse (at low stress) or an expanding crack-like rupture (at high stress). Because expanding pulses lead to an increasing net slip with increasing propagation distance, secondary nucleation is observed at the location of the perturbation at intermediate background stress.

3 General stability criterion

The numerical results clearly indicate that steady-state slip pulses driven by thermal pressurisation are unstable. How general is this result? In this Section, an approximate equation of motion for dynamic pulse-like ruptures is established, and utilised to determine a general stability criterion for steady-state slip pulses.

3.1 Slip pulse elastodynamics

In the elastodynamic equilibrium equation

[TABLE]

the stress redistribution functional has a form (Cochard and Rice, 1997)

[TABLE]

The kernel $M(u)$ assumes a simple form for anti-plane deformation:

[TABLE]

where $\mathcal{H}$ is the Heaviside function.

Consider a rupture in a form of a slip pulse of length $L(t)$ and total accumulated slip (dislocation) $b(t)$, the motion of which is specified by the coordinate of its tip, $x=\xi(t)$, advancing at generally non-uniform speed $v_{\mathrm{r}}(t)=\dot{\xi}(t)$. On spatial scales much larger than $L$, the pulse is seen as a singular dislocation

[TABLE]

where $\theta(x)$ is the arrival time of the pulse at the position $x$ (i.e., $\xi(\theta(x))=x$). Substituting this into (15) and (16), moving $\partial/\partial x$ under the integral in (16) yields after some manipulations (Appendix D)

[TABLE]

where

[TABLE]

with

[TABLE]

Terms $\phi_{\mathrm{Dirac}}$ and $\phi_{H}$ in the stress transfer functional $\phi=\phi_{\mathrm{Dirac}}+\phi_{H}$ correspond to the singular (Dirac) and the non-singular (step function $H$) terms in the slip gradient (Equation 19), respectively.

Since the region of applicability of “dislocation approximation” (20), $|x-\xi(t)|\gg L(t)$, excludes $x=\xi(t)$, the singular ($\delta_{\mathrm{Dirac}}$) term in (20) is of no consequence, and will be dropped in the foregoing, i.e.,

[TABLE]

3.2 Intermediate field of pulse

Let us introduce the coordinate frame moving with the dislocation (or the advancing front of the pulse),

[TABLE]

In the case of a steady (“steady state”) pulse motion, $\dot{b}=\dot{v}_{\mathrm{r}}=0$, one finds that $\bar{u}=v_{\mathrm{r}}/c_{\mathrm{s}}-X/c_{\mathrm{s}}(t-t^{\prime})$, and, in the case of anti-plane slip, (17), one recovers from (20) and (22) Weertman’s [1980] solution for a subsonic dislocation,

[TABLE]

In the general (non-steady) case, the Weertman’s solution (26) with instantaneous values of $b(t)$ and $v_{\mathrm{r}}(t)$ gives the leading order term in the near field of a dislocation (e.g. Eshelby, 1953; Markenscoff, 1980; Pellegrini, 2010). This near field can be defined by distances $X$ that are much smaller than a length scale $L_{\mathrm{out}}$ which characterises the unsteady motion of the pulse. For instance, if a dislocation accelerates or decelerates over a time scale $T_{\mathrm{out}}$ (e.g., $v_{\mathrm{r}}/\dot{v_{\mathrm{r}}}$ or $b/\dot{b}$), the associated length scale would be $L_{\mathrm{out}}=v_{\mathrm{r}}T_{\mathrm{out}}$. When considering a slip pulse, the approximation to a dislocation (Equation 20) only holds at distances much larger than the pulse length $L$, so that the near field of a dislocation corresponds to an intermediate field ($L\ll|X|\ll L_{\mathrm{out}}$) for a pulse, as long as the pulse “inner” lengthscale $L$ and the “outer” lengthscale $L_{\mathrm{out}}$ are separable,

[TABLE]

Furthermore, as shown by Eshelby (1953) on a particular example of accelerating dislocation motion, and by Markenscoff (1980); Callias and Markenscoff (1988); Ni and Markenscoff (2009) in the case of general motion $\xi(t)$ of a dislocation of invariant strength $b(t)=\text{constant}$, the next order term in the near field expansion of a non-uniformly moving dislocation is logarithmically singular. An extension of the results of Ni and Markenscoff (2009) for the $\phi_{\mathrm{Dirac}}$-expansion to the general case with arbitrary time-dependencies $\xi(t)$ and $b(t)$ yields (see Appendices E and F for details)

[TABLE]

Note that we nondimensionalised $X$ under the logarithm with “outer” lengthscale $L_{\mathrm{out}}$, but could have used other similar length for this purpose. This is due to the fact that any such scaling length contributes only to the high order terms, $O(X^{0})$, in the expansion. Admittedly, it would be advantageous to include these higher order terms to improve the approximation provided by this expansion (especially in view of the equation of motion discussed in the forthcoming). However, the actual expression for the $O(X^{0})$ correction is very cumbersome, and appears to depend on the history of slip (Callias and Markenscoff, 1988; Ni and Markenscoff, 2008), and, consequently, is not included in (28).

To find the near-field expansion for $\phi_{H}$ (Equation 23), we first write

[TABLE]

where

[TABLE]

with $\bar{u}$ defined in (24). Exploiting the similarity between the integrals in expressions for $\phi_{\mathrm{Dirac}}$ (Equation 22) and $\partial\phi_{H}/\partial t$ (Equation 30), and, in view of the $\phi_{\mathrm{Dirac}}$-expansion (equation 28), the leading term in the expansion for $\partial\phi_{H}/\partial t$ is given by

[TABLE]

Integrating, we have

[TABLE]

where, once again, the choice of normalising lengthscale under the logarithm ($L_{\mathrm{out}}$) is, apart from the order of magnitude considerations, somewhat arbitrary.

Combining (28) and (32), and simplifying, the near field expansion of the stress perturbation due to a moving dislocation takes the form

[TABLE]

3.3 Equation of motion of a moving dislocation

In view of (27), the stress $\tau(x,t)$ at intermediate distances from the pulse is approximately given by that of the steady-state dislocation with instantaneous strength $b(t)$, moving, as dictated by the steady-state pulse solution, at $v_{\mathrm{r}}\simeq v_{\mathrm{r,ss}}(b(t))$ within the “transient” background stress field, $\tau_{\mathrm{b,ss}}(b(t))$, i.e.,

[TABLE]

This type of approximation of the intermediate field of the unsteady dislocation appears to have been first suggested by Eshelby (1953, p. 251) when treating the particular example of a constant-strength dislocation accelerating from rest. Comparing (34) to (25) with (33), leads to an ordinary differential equation (ODE) describing the evolution of total slip $b(t)$ accrued in an unsteadily propagating pulse:

[TABLE]

where $v_{\mathrm{r}}=v_{\mathrm{r,ss}}(b(t))$ and $L=L_{\mathrm{ss}}(b(t))$ are the steady-state pulse velocity and width, respectively. In arriving to the form (35), the slowly space-varying $\ln|X/L_{\mathrm{out}}|$ term in (34) was approximated by its value at distances of few $L$ away from the trailing edge of the pulse. Equation (35) can be regarded as an “equation of motion” of an unsteady pulse, since once its solution $b=b(t)$ is known, the corresponding pulse trajectory follows by integration of $d\xi/dt=v_{\mathrm{r,ss}}(b(t))$.

In summary, the derived equation of motion relies on separation of spatial scales associated with the slip development within the pulse ($L$) and the evolution of the pulse net characteristics ($L_{\mathrm{out}}$), respectively, ($L\ll L_{\mathrm{out}}$). This scale separation allows to approximate unsteady pulse solution at a given instant of time by the steady-state solution (for a steadily propagating pulse) corresponding to the instantaneous value of total accrued slip $b(t)$, and other net pulse characteristics uniquely defined by the value of $b$ (i.e., $v_{\mathrm{r}}=v_{\mathrm{r,ss}}(b)$, $L=L_{\mathrm{ss}}(b)$, etc.). The evolution of the pulse “state variable” $b(t)$ is specified by the equation of motion (35).

3.4 Stability of steady-state pulse

Equation (35) allows to easily address the question of stability of a steady-state pulse solution (i.e., a solution of (35) with $db/dt=0$). If the rupture velocity of a steady-state pulse monotonically increases with total slip, $dv_{\mathrm{r,ss}}/db\geq 0$ and limited by $c_{\mathrm{s}}$ (see, for example, Garagash (2012) for steady rupture pulses driven by thermal pressurisation of pore fluid, and our Figure 3b,c), and in view of $L_{\mathrm{out}}\gg L$, the right hand side of (35) is a positive multiple of $db/dt$. It then follows from (35) that the sign of $db/dt$ is set by that of $\tau_{\mathrm{b}}-\tau_{\mathrm{b,ss}}(b)$, and, thus, a steady-state solution with $b(t)=b_{0}$ is stable to small perturbations if and only if the steady-state value of the background stress increases with slip, $(d\tau_{\mathrm{b,ss}}/db)_{|b=b_{0}}>0$. (Interestingly, a similar stability condition was cited by Rosakis (2001) without a proof).

For faults that dynamically weaken with slip, smaller levels of background stress are not inconsistent with larger required slip (and more pronounced weakening that comes with it) to drive a pulse rupture. We, therefore, expect the condition

[TABLE]

to be satisfied for a number of realistic constitutive laws (such as weakening by thermal pressurisation, as shown in Figure 7b), and thus inherently unstable steady-state pulse solutions. Indeed, initially steadily propagating slip pulses in a number of numerical studies utilising different models for the fault strength (e.g., Perrin et al., 1995; Beeler and Tullis, 1996; Zheng and Rice, 1998; Noda et al., 2009; Gabriel et al., 2012), eventually become unsteady, either growing (accelerating and accruing increasing levels of slip with distance travelled) or dying (shrinking and decelerating).

3.5 Perturbation growth rate

Let us write the equation of motion (35) in a shorthanded form $\tau_{\mathrm{b}}-\tau_{\mathrm{b,ss}}(b)=\mu\Psi(b)\,db/dx$, where

[TABLE]

and, as before, $v_{\mathrm{r}}=v_{\mathrm{r,ss}}(b)$ and $L=L_{\mathrm{ss}}(b)$. Nondimensional function $\Psi(b)$ is positive when, e.g., the steady-state rupture velocity is increasing with increasing net slip, as in the case of steady-state pulses driven by thermal pressurisation. Regardless of the sign of $\Psi(b)$, any small perturbation $\Delta b_{\mathrm{ini}}=(b-b_{0})_{\mathrm{ini}}$ from the steady-state pulse propagation with $b=b_{0}$, will initially evolve with the propagated distance $x$ as

[TABLE]

where $\Psi$ and $d\tau_{\mathrm{b,ss}}/db$ are evaluated at the baseline state $b=b_{0}$. The exponential form (38) is qualitatively consistent with the numerical simulation using thermal pressurisation as a weakening mechanism (Figures 6 and 7a).

In the event when the slip perturbation is seeded by a background stress perturbation $\Delta\tau_{\mathrm{b}}$ localised in space over the dimension $\Delta x$, as is the case in our numerical simulations, the corresponding level of equivalent “initial” slip perturbation $\Delta b_{\mathrm{ini}}$ (that will persist and evolve according to (38) for $x>\Delta x$) can be estimated as

[TABLE]

The perturbation exponential growth rate, given by

[TABLE]

is therefore expected to be independent of the (small) perturbation amplitude. These general observations are consistent with the numerical simulations, which show a linear scaling between the perturbation in stress and the resulting slip perturbation, and the independence of the growth rate on the stress perturbation amplitude (Appendix C, Figure 14).

The steady-state solutions presented in Section 2 provide the relationships between $\tau_{\mathrm{b,ss}}$, $v_{\mathrm{r,ss}}$ and slip $b$ (see Figure 3) required to compute a theoretical estimate of the growth rate using Equation (40), leaving only the ratio $L_{\mathrm{out}}/L$ as an unconstrained parameter. Using $L_{\mathrm{out}}/L=10$ or a constant $L_{\mathrm{out}}/L^{*}=10$ and the steady-state pulse width $L=L_{\mathrm{ss}}$ produces the results shown in Figure 11 (dotted and solid lines, respectively), where the growth rates estimated from numerical simulations are also displayed for comparison. The agreement between theoretical and numerical estimates with either choice for $L_{\mathrm{out}}$ is very satisfactory, and illustrates the applicability of the pulse equation of motion (35). Since the ratio $L_{\mathrm{out}}/L$ only appears in the logarithmic term, the resulting growth rate is not very sensitive to the specific choice for this unconstrained quantity. It appears that choosing $L_{\mathrm{out}}$ to be several times larger than $L$ produces reasonable predictions, consistent with the assumption (27).

3.6 Validity of equation of motion

The key underlying assumption in our derivation of the approximate equation of motion for the slip pulse (Equation 35) is that the pulse is in “quasi-steady-state”, i.e., its characteristics (accrued slip $b$, length $L$ and speed $v_{\mathrm{r}}$) change slowly on the timescale of slip. This assumption can be translated in terms of propagation distance, since $b,L$ and $v_{\mathrm{r}}$ do not vary appreciably over propagation distances of the order of the pulse length $L$: the quasi-steady-state approximation is then valid as long as $|d(b/\delta_{\mathrm{c}})/d(x/L^{*})|\ll 1$, that is, $(\mu/\tau_{0})|db/dx|\ll 1$. This assumption can be validated from the equation of motion itself: indeed, using the notation introduced in Equation (37),

[TABLE]

which is a function of slip $b$ and stress $\tau_{\mathrm{b}}$, plotted in Figure 12. For elevated background stresses, around $\tau_{\mathrm{b}}/\tau_{0}=0.9$, the normalised slip gradient $(\mu/\tau_{0})(db/dx)$ remains significantly less than unity. This is also the case throughout the regime of growing pulses (i.e., when $\tau_{\mathrm{b}}>\tau_{\mathrm{b,ss}}$). We therefore expect the equation of motion to provide an adequate description of the pulse dynamics under those conditions. For arresting pulses, the assumption of quasi-steady-state becomes invalid as slip decreases, with the magnitude of the slip gradient rapidly becoming of the order of unity, notably under low background stresses. This can be understood by considering that steady-state pulses associated with small slip correspond to elevated background stresses and low rupture speeds (Figure 3b,c): the regime of arresting pulses under very low stresses $\tau_{\mathrm{b}}\ll\tau_{\mathrm{b,ss}}$ is therefore too far from steady-state and the pulse is expected to arrest quickly compared to the duration of slip.

4 Discussion and implications

The pulse equation of motion and stability analysis demonstrate that steady-state pulses are unstable if, e.g., $d\tau_{\mathrm{b,ss}}/db<0$ and $dv_{\mathrm{r}}/db>0$, or more generally if the exponential growth rate is positive, which, in view of (40) and (37), corresponds to

[TABLE]

This condition is satisfied for pulses driven by thermal pressurisation, and the numerical simulations confirm qualitatively and quantitatively this instability. In the light of these results, two key questions arise: what do steady-state pulse solutions tell us about the dynamics of rupture in general? What does the existence of unstable slip pulses imply for earthquake dynamics and strength of faults?

4.1 Significance of steady-state pulse solutions

Steady-state slip pulses arise spontaneously in fully dynamic rupture simulations when the nucleation and background stress conditions are at the transition between arresting and growing ruptures (Noda et al., 2009; Schmitt et al., 2011; Gabriel et al., 2012). Therefore, the conditions leading to the existence of steady-state solutions coincide with those allowing for the existence of sustained ruptures. In other words, steady-state pulse solutions inform us about the overall “strength” of an interface, in the sense that they provide us with the critical conditions required for ruptures to propagate beyond their nucleation patch.

Our results complement the framework provided by Zheng and Rice (1998) who determined the critical background stress level ($\tau_{\mathrm{pulse}}$) separating the regime of exclusively pulse-like ruptures under low stress conditions and the regime where both crack and pulse rupture modes are possible under high stress conditions. Here we show both theoretically and numerically (for the case of thermal pressurisation) that the pulse mode of rupture exists within the entire range of background stress (low and high), while the dynamics of the pulse (spontaneous decay leading to arrest or spontaneous growth leading to either transition into crack-like rupture or nucleation of a secondary rupture in the pulse wake) can be extracted from a steady-state pulse analysis, like that conducted by Garagash (2012) and summarised in Section 2.1 for the case of thermal pressurisation. Although we do not establish conditions for prevalence of the pulse-like mode for ruptures driven by thermal pressurisation (such as the $\tau_{\mathrm{pulse}}$ threshold of Zheng and Rice (1998) for velocity-weakening friction case), we suspect that the $\tau_{\mathrm{pulse}}$ threshold in this case would correspond to the minimum level of background stress at which the secondary rupture nucleated in the wake of growing primary pulse is crack-like. Therefore, the solution to the steady-state pulse problem associated with a particular constitutive behaviour provides a tool to determine the exact conditions (notably in terms of background stress) leading to the existence of sustained ruptures. Our analysis on the role of unstable slip pulses in controlling the growth of large scale ruptures is consistent with recent theoretical results from Brener et al. (2018), who analysed numerically the stability of slip pulses driven by a nonlinear rate-dependent friction law. Numerical simulations indicate that such slip pulses are also unstable to small perturbations, and Brener et al. (2018) argue that such instabilities can be viewed as the nucleation process of large ruptures.

In the case of thermal pressurisation, the minimum dynamic strength is zero and thus there is no lower stress limit for the existence of dynamic steady-state slip pulses. Therefore, faults governed by thermal pressurisation have theoretically no “strength”: thermal pressurisation allows for large enough pulses to propagate regardless of the initial background stress. However, theoretical slip pulses propagating under very low stress conditions bear large slip and slip rate, and require nucleation conditions characterised by either very high local stresses or large nucleation region (with modestly elevated stress). The question of the minimum stress required for ruptures to grow is therefore linked to the nucleation conditions of those ruptures. This was illustrated by Gabriel et al. (2012) in the context of a slip rate dependent constitutive law, who showed that the threshold background stress between arresting and growing pulses (i.e., the steady-state pulse regime) scales with the size of the nucleation patch used in their simulations. The nucleation conditions probably enforce the selection of a specific characteristic pulse width, stress drop and slip rate, and the background stress level outside the nucleation patch selects whether the rupture will become crack-like or an expanding pulse or decaying pulse, the boundary between the latter two regimes being determined by the steady-state stress for that pulse.

4.2 Complexity of earthquake ruptures

Our results show that unstable slip pulses can produce remarkably complex rupture events, even when the background stress and conditions are uniform (except for an infinitesimal perturbation). Rupture complexity is thus not systematically linked to complexity or heterogeneity in stress or strength conditions, but arise spontaneously when ruptures propagate as slip pulses.

A notable feature observed for pulses arresting due to negative stress perturbations is that the arrest is abrupt. Negative perturbations in slip, slip rate and pulse width grow exponentially over time until the rupture stops. This is best illustrated by computing the (one dimensional) moment rate,

[TABLE]

which is shown in Figure 13 for a pulse propagating at $\tau_{\mathrm{b}}/\tau_{0}=0.7$ and arresting due to a small negative perturbation. The moment rate is initially constant, corresponding to the steadily propagating pulse solution. Upon encountering the stress perturbation, moment rate decreases exponentially and drops abruptly to zero. Such rapid variations in moment rate are responsible for the radiation of high frequency waves in the far field, and we confirm here that such high frequencies associated with sudden rupture arrest can arise without strong stress or strength heterogeneities on the fault plane. Similar conclusios were established by Cochard and Madariaga (1994) and Gabriel et al. (2012) in simulations using strong velocity weakening friction.

Another key feature associated with the slip pulse instability is the transition from pulse-like to crack-like rupture due to positive stress perturbations under high uniform background stress (Figures 8 and 9). Here again, we observe a remarkable complexity that emerges spontaneously in the absence of any preexisting fault heterogeneities. Above a critical background stress (here $\tau_{\mathrm{b}}/\tau_{0}\gtrsim 0.79$), the slip pulse transitions directly to an expanding crack. Under low stress conditions, numerical results indicate that growing ruptures become expanding pulses. One important consequence of the transition to expanding pulse is that as the pulse further propagates, the accrued slip grows approximately linearly with propagation distance and therefore we expect a logarithmic stress build-up near the starting point of growth. This is observed in our simulations when the background shear stress is close to the threshold for the direct transition into crack-like rupture (see Figures 9, 10).

What this transition illustrate is that unstable steady-state pulses evolve towards the most stable rupture mode, either self-similar pulse or expanding crack, according to the current background stress level. However, a peculiarity exhibited by our results is that rupture arrest is also a strong attractor (when perturbations are negative), so that a nascent slip pulse propagating in an overall high stress regime might arrest on its own if negative stress perturbations are encountered.

5 Conclusions

We performed numerical simulations and a theoretical analysis that demonstrate that steady-state slip pulses are unstable if the accrued slip (“dislocation”) is a decreasing function of the background stress, i.e., $d\tau_{\mathrm{b}}/db\leq 0$. This instability condition is satisfied for slip pulses driven by thermal pressurisation of pore fluids. During instability, slip, slip rate and pulse width perturbations grow exponentially. If the initial stress perturbation leading to instability is negative, ruptures eventually arrest in an abrupt manner; conversely, if the stress perturbation is positive, rupture mode changes and transitions to a growing pulse (at low stress) or an expanding crack (at high stress). The growth rate of perturbations is predicted quantitatively by an approximate equation of motion for a dislocation with variable net slip (Equation 35).

The regime of steady-state pulse solutions appears naturally in dynamic rupture simulations at the transition between spontaneously expanding ruptures (growing pulses) and spontaneously arresting ruptures, at a stress level that depends on the nucleation conditions. Once nucleation conditions are established, the steady-state pulse solution provides a stress limit below which ruptures will spontaneously stop, which is best considered as the “strength” of the interface (e.g., Lapusta and Rice, 2003; Rubinstein et al., 2004; Noda et al., 2009).

The unstable character of steady-state slip pulses generates a remarkable complexity of ruptures, including abrupt arrest, pulse to crack transitions and secondary rupture nucleation in the wake of a propagating pulse, even though stress and material parameters are homogeneous (nonwithstanding an infinitesimal perturbation) along the fault. Pulse-like ruptures seem to be the main rupture mode for many crustal faults (Heaton, 1990), and it is therefore expected that earthquake dynamics along these faults is driven at least in part by spontaneous instabilities. One key consequence is that abrupt arrest of ruptures may not be the signature of strong preexisting stress of strength heterogeneities along faults. At this stage, it remains to be explored how slip pulse instabilities evolve along heterogeneous faults, and further work in this direction is currently conducted. Preliminary simulations suggest that transient pulses (i.e., non steady-state) continuously grow or shrink as they cross regions of high and low stress, respectively, their eventual arrest being dictated by finite amplitude stress perturbations.

Acknowledgements.

All authors contributed equally to the work, and are listed in alphabetical order. This work was supported by the UK Natural Environment Research Council through grant NE/K009656/1 to NB, by the Canada Natural Science and Engineering Research Council though Discovery grant 05743 to DIG and by MEXT KAKENHI Grant Number 26109007 to HN. The results of this paper can be reproduced by direct implementation of the analytical formulae and numerical methods described in the main text and appendices.

Appendix A Convolution kernels for fault strength and temperature

The convolution kernel $\mathcal{K}$ is given by (Rice, 2006; Garagash, 2012)

[TABLE]

if $\chi\neq 1$, or by the limit of that expression as $\chi\rightarrow 1$ if $\chi=1$. The function $\mathcal{A}$ depends on the spatial distribution of strain rate across the fault, and for our choice of a Gaussian distribution we have

[TABLE]

The pore pressure and temperature evolution on the fault plane ($y=0$) can be computed from the strength evolution as

[TABLE]

and

[TABLE]

Appendix B Numerical methods

B.1 Steady-state problem

The method of solution for the steady-state pulse is the same as that employed by Platt et al. (2015) in a more complex case (including thermal decomposition in addition to thermal pressurisation), and reviewed by Viesca and Garagash (2018). For completeness, we provide here a description of the technique in the simple case of a finite pulse driven by thermal pressurisation only.

Normalising the slip by $\delta_{\mathrm{c}}$, time by $T^{*}$, stresses by $\tau_{0}$, distances by $\bar{\mu}\delta_{\mathrm{c}}/\tau_{0}$ and slip rate by $\delta_{\mathrm{c}}/T$, we rewrite the governing Equations (1) and (8) as

[TABLE]

and

[TABLE]

where normalised variables are denoted by a tilde, $\mathcal{H}$ is the Heaviside function, and $\chi=\alpha_{\mathrm{hy}}/\alpha_{\mathrm{th}}$. In Equation (49), we changed the integration variable from time to space by noting that $\tilde{t}=\tilde{x}\tilde{T}/\tilde{L}$. The integrals in (48) and (49) are further normalised using the transformed space coordinate $y=2\tilde{x}/\tilde{L}-1$, which results in

[TABLE]

The condition (3) is similarly rewritten as

[TABLE]

The idea now is to approximate the above integrals with Gauss-Chebyshev quadratures (Viesca and Garagash, 2018). Because we expect the slip rate $\tilde{V}(y)$ to behave as $\sqrt{1\pm y}$ near $y\mp 1$ (i.e., square-root behaviour of the slip profile near the rupture tip and tail), we introduce the function $v(y)$ as

[TABLE]

which becomes the unknown (regular) function we are looking to approximate. Using the approximations

[TABLE]

with

[TABLE]

for $i=1,\ldots,n$, $j=1,\ldots,n+1$, and

[TABLE]

with

[TABLE]

for $p=1,\ldots,n$, the governing equations become a linear system:

[TABLE]

where $\tilde{\tau}_{i}=\tilde{\tau}(Y_{i})$, $\tilde{\tau}_{j}=\tilde{\tau}(y_{j})$ and $v_{j}=v(y_{j})$. In the system above, the normalised stress $\tilde{\tau}$ needs to be differentiated with respect to $y$, and evaluated at both sets of points $y_{j}$ and $Y_{i}$. Given the knowledge of the set of $\tilde{\tau_{i}}$, we use barycentric interpolation and Chebyshev differentiation matrices to compute

[TABLE]

where $L_{ji}$ is an interpolation matrix (Viesca and Garagash, 2018) and $D_{pj}$ is a Chebyshev differentiation matrix (Trefethen, 2000), and we sum over repeated indices. In summary, we arrive at the following linear system:

[TABLE]

where $K_{ij}=w_{j}/(\pi(Y_{i}-y_{j}))$ and $S_{ij}=w_{j}\mathcal{H}(Y_{i}-y_{j})\mathcal{K}((Y_{i}-y_{j})T/2;\chi)$. Equating (61) and (62), we obtain a total of $n+2$ equations, with $n+2$ unknowns that are $v_{j}$ ($j=1,\ldots,n$), $\tilde{L}$ and $\tilde{T}$. This system is solved using the Newton-Raphson iterative algorithm.

B.2 Expression of stress transfer functional

Consider a spatial domain of length $\lambda$. Let $D_{p}(t)$ and $\dot{D}_{p}(t)$ denote the spatial discrete Fourier transform coefficients of the slip and slip rate perturbations, respectively, where indices $p$ correspond to wavenumbers $k_{p}=2\pi p/\lambda$. The discrete Fourier transform coefficients of the stress transfer functional $\phi$ are given by (Perrin et al., 1995)

[TABLE]

where $W(u)=\int_{0}^{\infty}\left(J_{1}(x)/x\right)dx$, and $J_{1}$ is the Bessel function of the first kind of order one. An inverse Fourier transform of $F_{p}(t)$ provides the value of $\phi$ in the space-time domain.

B.3 Dynamic problem

The technique employed to solve the elastodynamic problem is essentially following the spectral boundary integral method of Lapusta et al. (2000), adapted to our specific choice of constitutive behaviour (thermal pressurisation with constant friction coefficient). In this method, the dynamic stress transfer functional is evaluated in the Fourier domain, taking advantage of the efficiency of Fast Fourier Transform (FFT) algorithm.

The space domain is discretised into nodes $x_{i}=ih$, $i=1,\ldots,N$. Time is discretised into steps $t_{n}$, $n=0,\ldots,N_{t}$, with a constant spacing $\Delta t$. We denote with subscripts $i$ and superscripts $n$ the discretised variables at node $(x_{i},t_{n})$.

We first determine a steady-state solution for a uniform background stress and a given diffusivity ratio. The stress and slip rate distributions, $\tau_{\mathrm{ss}}$, $V_{\mathrm{ss}}$, are interpolated onto our regular grid at each node $(x_{i},t_{n})$, so that $\tau_{\mathrm{ss},i}(t_{n})$ and $V_{\mathrm{ss},i}(t_{n})$ are precomputed and stored a priori. At time $t_{0}$, we initialise the perturbations in slip ($\Delta\delta_{i}$), slip rate ($\Delta V_{i}$), stress ($\Delta\tau_{i}$) and strength ($\Delta\tau_{\mathrm{f},i}$) with zeros at all nodes.

Let us consider that all variables are known at a given time step $t_{n}$, including the entire slip rate perturbation history (and its Fourier coefficients, for use in the spectral boundary integral algorithm). The computation of variables at time step $t_{n+1}=t_{n}+\Delta t$ is conducted as follows:

1

Make a first estimate of the slip perturbation assuming a slip rate perturbation equal to that at time step $t_{n}$:

[TABLE]

2

Estimate the perturbation in dissipation rate (denoted $\Delta(\tau V)$) for the interval $[t_{n},t_{n+1}]$ as

[TABLE]

and compute the perturbation in strength as

[TABLE]

which corresponds to a mid-point approximation of the integral in (14). The computation of (67) requires the storage of the full history in $\Delta(\tau V)$.

3

Compute the Fourier coefficients $D_{p}^{*}$ and $\dot{D}_{p}^{*}$ of the first estimates of slip and slip rate perturbation profiles at time step $t_{n+1}$, where subscripts $p$ indicate wavenumber indices. This operation is performed using the FFT algorithm. Then estimate the stress transfer functional in the Fourier domain as (see Equation (64))

[TABLE]

where $W_{p}^{k}=W(|k_{p}|c_{\mathrm{s}}t_{k})$ (see Appendix B.2). Using an inverse FFT, compute an estimate $\phi^{*}_{i}$ of the stress transfer functional.

4

Compute the total strength $\tau^{*}_{\mathrm{f},i}=\tau_{\mathrm{f,ss},i}+\Delta\tau^{*}_{\mathrm{f},i}$, and the total stress $\tau_{\mathrm{stuck},i}$ that would be applied if slip rate was zero, given by

[TABLE]

Slip rate is nonzero where $\tau^{*}_{\mathrm{f},i}>\tau_{\mathrm{stuck},i}$. At those nodes, assign $\Delta\tau^{*}_{i}=\Delta\tau_{\mathrm{f},i}^{*}$, and compute the slip rate perturbation as

[TABLE]

Where $\tau^{*}_{\mathrm{f},i}<\tau_{\mathrm{stuck},i}$, assign $\Delta\tau^{*}_{i}=\tau_{\mathrm{stuck},i}-\tau_{\mathrm{ss},i}$ and $\Delta V_{i}^{*}=-V_{\mathrm{ss},i}^{n+1}$.

5

Repeat steps 1 to 4 using $(\Delta V_{i}^{*}+\Delta V_{i}^{n})/2$ and $(\Delta\tau_{\mathrm{f},i}^{*}+\Delta\tau_{\mathrm{f},i}^{n})/2$ instead of $\Delta V_{i}^{n}$ and $\Delta\tau_{\mathrm{f},i}^{n}$, respectively. The convolutions in Equations (67) and (68), which are the most computationally intensive steps, are not recomputed entirely but simply updated because only the last term has changed. The resulting slip, slip rate, stress and strength perturbations are the final predictions at the next time step $\Delta\delta_{i}^{n+1}$, $\Delta V_{i}^{n+1}$, $\Delta\tau_{\mathrm{f},i}^{n+1}$ and $\Delta\tau_{i}^{n+1}$, respectively.

Appendix C Perturbation amplitude

Figure 14 shows the time evolution of slip perturbations following negative perturbations in background stress of $10^{-4}$, $10^{-3}$ and $10^{-2}$ in amplitude. In all simulations the reference background stress is $\tau_{\mathrm{b}}/\tau_{0}=0.7$ and the diffusivity ratio is $\alpha_{\mathrm{hy}}/\alpha_{\mathrm{th}}=1$. The growth of the slip perturbation is exponential, and the growth rate does not depend on the amplitude of the stress perturbation. The initial jump in normalised slip (and other normalised variables, see Figure 6) is directly proportional to the amplitude of the stress perturbation.

Appendix D Expression for $\phi_{H}(x,t)$

The contribution $\phi_{H}$ to the stress-transfer functional $\phi$ (Equation 16) from the second term in the expression for the slip gradient (Equation 19) can be written, after moving $\partial/\partial x$ under the integral and substituting $dx^{\prime}=(d\theta(x^{\prime})/dx^{\prime})^{-1}d\theta$, in the following form:

[TABLE]

where $M^{\prime}(u)=dM/du$. Changing the order of integration in the above double integral

[TABLE]

and carrying out the integral in $t^{\prime}$, one finds a single-integral expression for $\phi_{H}(x,t)$. This expression, after changing the integration variable symbol from $\theta$ for $t^{\prime}$, is given in the main text (Equation 23 with 24).

Appendix E Deduction of expression (28) for $\phi_{\mathrm{Dirac}}(x,t)$

based on the work of Ni and Markenscoff [2009]

Equation (28) can be established by accounting for the time-dependence of slip $b(t)$ in the derivation of the results of Ni and Markenscoff (2009) (their Equation (5.18)), who only considered dislocations with constant slip $b$. In practice, our Equation (28) results from carrying out the time-derivative of $b(t)$ from Ni and Markenscoff’s Equation (3.15) to obtain a more general form of their Equation (3.16), and then equating their Equation (5.16) to the modified Equation (3.16).

It appears that Ni and Markenscoff (2009), and other references of Markenscoff and co-workers give different sign (minus) in front of the logarithmic term compared to the result used here in (28). A negative sign in front of the logarithmic term is inconsistent with Eshelby’s 1953 example (uniform acceleration), although one should be cautioned that there is a typographical error in the definition of $s_{0}$ used to evaluate the stress expansion (his Equation (15)) in Eshelby’s paper. It should read $s_{0}=|x_{0}-\xi(t)|/c$ instead of $s_{0}=\sqrt{t^{2}-[x-\xi(t)]^{2}}$ given directly under his Equation (15). This definition of $s_{0}$ has to be corrected in order to evaluate the logarithmic term in the stress-expansion correctly. Markenscoff’s negative sign in front of the logarithmic term is also inconsistent with examples of direct numerical evaluation of $\phi_{\mathrm{Dirac}}$ for prescribed dislocation evolution functions $\xi(t)$ and $b(t)$ (computations performed using Mathematica). Furthermore, we explicitly derive the $\ln$-term in the case of constant $v_{\mathrm{r}}$ and $b=b(t)$ below (Appendix F), which supports the sign used here. On the basis of the arguments above, we conclude that there is a typographical sign error in the work of Markenscoff and co-workers.

Appendix F Alternative derivation of $\phi_{\mathrm{Dirac}}(x,t)$ for particular case

of a steady motion of dislocation of variable strength $b=b(t)$

Introducing the coordinate $X$ moving with the crack tip $X=\xi(t)-x$, and using the following notation $\Delta t=t-t^{\prime}$ and $\Delta\xi=\xi(t)-\xi(t^{\prime})$, we rewrite $\phi_{\mathrm{Dirac}}$ as

[TABLE]

with

[TABLE]

The integral in (73) can be decomposed in that over $\Delta t$ within and outside a $t_{\epsilon}$-window. In the former, the dynamic quantities (dislocation and rupture velocity) can be approximated using their current rates, i.e.,

[TABLE]

The part of the integral (73) for $\Delta t>t_{\epsilon}$ is bounded by $O(b/t_{\epsilon})$.

Thus, we focus on the integral for $\Delta t<t_{\epsilon}$ with expectation that it provides the singular part of the near-field ($X\rightarrow 0$) of $\phi_{Dirac}$. Furthermore, let us restrict the consideration to steady dislocation motion $\ddot{\xi}=0$, which simplifies the variable $\bar{u}$ dependence on $\Delta t$ to the following

[TABLE]

where we have renamed $\dot{\xi}=v_{\mathrm{r}}$)

With the above, and changing integration variable to $\bar{u}$, i.e., $\Delta t=\Delta t(\bar{u})$, we have

[TABLE]

and the “lower” bound of integration corresponds to minimum value of $\Delta t$ given by $\Delta t_{1}=-X/(c_{s}-\dot{\xi})$ for which the integrand is non-zero, i.e., $\bar{u}(\Delta t_{1})=1$, while the “upper” bound $\bar{u}_{\epsilon}=\bar{u}(t_{\epsilon})$. Explicit integration, expanding in series in small $X$ and retaining the singular ($1/X$ and $\ln|X|$) terms leads to

[TABLE]

The first in the above corresponds to the steady dislocation with constant strength, while the second is the correction for variable strength.