The stability of frictional sliding on dip-slip and finite length faults

SUMMARY

This paper examines the linear stability of sliding on faults embedded in a 2-D elastic medium that obey rate and state friction and have a finite length and/or are near a traction-free surface. Results are obtained using a numerical technique that allows for analysis of systems with geometrical complexity and heterogeneous material properties; however only systems with homogeneous frictional and material properties are examined. Some analytical results are also obtained for the special case of a fault that is parallel to a traction-free surface. For velocity-weakening faults with finite length, there is a critical fault length |$L^{*}$| for unstable sliding that is analogous to the critical wavelength |$h^{*}$| that is usually derived from infinite fault systems. Faults longer than |$L^{*}$| are linearly unstable to perturbations of any length. On vertical strike-slip faults or faults in a full-space |$L^{*} \approx h^{*}/e$|⁠, where e is Euler’s number. For dip-slip faults near a traction-free surface |$L^{*} \le h^{*}/e$| and is a function of dip angle |$\beta$|⁠, burial depth d of the fault’s up-dip edge and friction coefficient. In particular, |$L^{*}$| is at least an order of magnitude smaller than |$h^{*}$| on shallowly dipping (⁠|$\beta < 10^\circ$|⁠) faults that intersect the traction-free surface. Additionally, |$L^{*} \approx h^{*}/e$| on dip-slip faults with burial depths |$d \ge h^{*}$|⁠. For sliding systems that can be treated as a thin layer, such as landslides, glaciers or ice streams, |$L^{*} = h^{*}/2$|⁠. Finally, conditions are established for unstable sliding on infinitely-long, velocity-strengthening faults that are parallel to a traction-free surface.

1 INTRODUCTION

How does the geometry of fault systems affect the stability of frictional sliding? Most knowledge of frictional stability comes from analyses of spring-slider systems (e.g. Ruina 1983; Gu et al. 1984; Ranjith & Rice 1999) or faults embedded in an infinite elastic medium (e.g. Rice & Ruina 1983; Dieterich 1992; Rice et al. 2001; Uenishi & Rice 2003; Rubin & Ampuero 2005; Ampuero & Rubin 2008). These two systems do not include important aspects of fault geometry such as proximity to a traction-free surface, fault dip angle or finite fault length in many cases. This paper examines the linear stability of sliding on finite length faults that obey rate and state friction (RSF) and are embedded in a 2-D elastic continuum with homogeneous material properties. As in many previous studies, the focus here is on the quasi-static regime where inertia is neglected (e.g. Rice & Ruina 1983; Uenishi & Rice 2003; Viesca 2016a, b; Aldam et al. 2017; Heimisson et al. 2019; Ozawa et al. 2024).

The specific finite length geometries considered are: a fault in an infinite full-space; a fault parallel to a traction-free surface and a dip-slip fault or vertical strike-slip fault in a half-space. Including these features provides more accurate assessments of the sliding stability of natural fault systems. First, because all faults have a finite length and secondly, because many faults are near the surface of the earth or the seafloor. Additionally, landslides (Handwerger et al. 2016), ice streams (Lipovsky & Dunham 2017) and glaciers (Zoet et al. 2020) also exhibit sliding behaviour that can be described by frictional processes. The results show that these basic geometrical effects cause significant departures from long-standing results on linear stability behaviour.

For a fault in an elastic medium, the stability of sliding can be assessed by considering a balance between the rate at which elastic stress stored in the fault system can be unloaded, and the rate at which shear stress on the fault is reduced (i.e. fault weakening) in response to sliding (e.g. Scholz 2019). Unstable sliding initiates when the fault weakening rate is higher than the elastic unloading rate. For frictional sliding, the changes in shear stress on a fault are described by the RSF equations wherein the evolution of the friction coefficient |$\mu$| on a surface is a function of the sliding rate v and an internal state variable |$\theta$|

where |$\mu _0$| is a reference friction coefficient for steady state sliding at the reference velocity |$v_0$|⁠, a and b are constitutive parameters, and |$d_c$| is a critical slip distance that is related to the amount of slip needed to attain a new steady state after changes in sliding velocity (Dieterich 1979; Ruina 1983; Marone 1998). The reference values |$\mu _0$| and |$v_0$| are not arbitrary, but are determined from experimental measurements of how the steady state friction coefficient of a given material depends on the sliding velocity. However, pairs of values (⁠|$\mu _0$|⁠, |$v_0$|⁠) can be chosen from among a set of such measurements. The evolution of the state variable is commonly described using the aging or slip laws:

For velocity-weakening friction (⁠|$a-b < 0$|⁠; or |$a/b < 1$|⁠), fault weakening will occur due to a reduction in the friction coefficient as sliding rate increases. On deformable faults, weakening can also occur due to a coupling between fault slip and changes in normal stress, which can lead to unstable behaviour for velocity-strengthening friction below some critical value of |$a/b$|⁠. This effect has been shown to exist on bimaterial and poroelastic faults (Rice et al. 2001; Ranjith 2014; Heimisson et al. 2019), on faults with fault-valve behaviour (Ozawa et al. 2024), and on faults that lack geometric reflection symmetry across the sliding surface (Aldam et al. 2016). Lack of geometric reflection symmetry is a very general feature of fault systems. So too then is the possibility of unstable behaviour on velocity-strengthening faults. The results in this paper take a first step in establishing the range of parameters where this behaviour occurs on faults that are near a traction-free surface.

Before proceeding, it will be helpful to differentiate between linear and nonlinear methods of stability analysis, as well as establish some terminology related to frictional sliding stability. Studies of frictional sliding stability can be grouped according to the use of linear analysis methods (e.g. Rice & Ruina 1983; Ruina 1983; Rice et al. 2001; Heimisson et al. 2019; Ozawa et al. 2024) or nonlinear analysis methods (e.g. Gu et al. 1984; Dieterich 1992; Ranjith & Rice 1999; Uenishi & Rice 2003; Rubin & Ampuero 2005; Ampuero & Rubin 2008; Viesca 2016a; Ray & Viesca 2017). In a linear analysis (the subject of this paper), boundary conditions are defined such that a spatially uniform steady state exists on the fault and the equations governing fault slip are simplified through linearization about this steady state. In addition, the initial conditions must be defined to reside within the neighbourhood of the steady state in order for the linearized governing equations to remain valid. A stability analysis of the linearized equations determines if a set of initial conditions, that includes small perturbations to the steady state sliding velocity, will grow or decay. If the fault system is linearly stable then the sliding velocity will return to the steady state after it is perturbed. If the fault system is linearly unstable then the sliding velocity will grow exponentially until the linearized equations are no longer valid, at which point the nonlinear governing equations must be considered to determine the subsequent evolution. Throughout this paper the words ‘instability’ and ‘unstable’ are used according to this meaning. The word ‘perturbation’ is also used frequently and will always refer to a deviation from steady state.

Nonlinear stability analysis concerns the full nonlinear governing equations. Various initial and boundary conditions can be considered since there are no restrictions imposed as in a linear analysis. In a nonlinear analysis, ‘instability’ usually means that slip speeds would become infinitely large if they were not restricted by inertia (e.g. Rubin & Ampuero 2005; Ampuero & Rubin 2008). Because the results in this paper only concern behaviour in the linear regime, the nonlinear meaning of ‘instability’ is not used.

Linear stability analysis of spring-slider systems leads to the classical concept of a critical patch length, that is interpreted as the minimum length of a slipping fault patch that is required for an unstable sliding event to develop (e.g. Dieterich ; Rice 1993; Scholz 1998). The critical patch length |$l_c$| is obtained by equating the critical stiffness of an RSF spring-slider system to the stiffness of a crack embedded in an elastic medium, such that

where |$G^{\prime } = G$| for antiplane sliding and |$G^{\prime } = G/(1 - \nu )$| for in-plane sliding, with G the shear modulus and |$\nu$| the Poisson ratio. In eq. (3), |$\eta$| is a numerical pre-factor of order unity that depends on the assumed stress and strain conditions for computing the stiffness of the elastic crack (see table 1 in Dieterich 1992).

However, the stiffness of a fault in an elastic continuum is a quantity that evolves through space and time as the fault slips (e.g. Rice & Ruina 1983; Horowitz & Ruina 1989), whereas eq. (3) is obtained by assuming the existence of a unique critical stiffness. Linear stability analysis of infinitely long faults in elastic continua leads to a critical wavelength, rather than a critical stiffness or patch size. The critical wavelength |$h^{*}$| can be found analytically for the special case of an infinitely long fault with constant frictional properties and effective normal stress |$\sigma$|⁠, embedded in an infinite, 2-D elastic full-space with homogeneous properties. Allowing the fault to be infinitely long simplifies the mathematical analysis sufficiently to obtain an equation for |$h^{*}$| (Rice & Ruina 1983; Rice et al. 2001):

Here the symbol |$h_F^{*}$| is used to denote the special value of |$h^{*}$| for a homogeneous fault in a full-space. Because the derivation of |$h^{*}_F$| involves an infinitely long fault, its proper definition is the critical wavelength of an infinitely long, sinusoidal perturbation to the steady state.

The fact that the critical wavelength |$h^{*}_F$| is identical to the critical patch length |$l_c$| (to within the numerical constant |$\eta$|⁠) naturally leads to the conclusion that unstable sliding on velocity-weakening, elastically deformable faults can only occur when perturbed slipping patches are longer than |$l_c$|⁠. However, this is not true for sliding behaviour in the linear regime, because any perturbation that contains wavelengths larger than |$h^{*}_F$| will grow exponentially until the linearized equations are no longer valid. Importantly, since an initial perturbation need not have a sinusoidal form, this does not impose any conditions on the spatial extent of the perturbation (see Section 2.1.1).

When faults are not infinitely long or when a spring-slider model cannot capture important features of a fault system, then it becomes difficult to apply analytical methods of linear stability analysis. In this paper these difficulties are overcome by using a numerical method for conducting linear stability analysis of 2-D finite length fault systems. The method can be applied to any fault system for which stress change functions are available (defined in the next section), and can accommodate features such as heterogeneous material properties or multiple faults.

The results of this paper show that the linear stability behaviour of finite length faults is properly characterized by a critical fault length |$L^{*}$|⁠, the value of which is related to the critical wavelength of a corresponding infinitely long fault. Once the fault length is larger than |$L^{*}$|⁠, then a perturbation of any size will grow exponentially within the linear regime. Throughout this paper, values of |$h^{*}$| are referred to as ‘critical wavelengths’ and the symbol |$L^{*}$| is used to denote a ‘critical fault length’. Subscripts are used for both |$h^{*}$| and |$L^{*}$| to differentiate between specific geometries. On vertical strike-slip faults or faults in a full-space, |$L_F^{*} \approx h^{*}_F/e$| where e is Euler’s number. For dip-slip faults near a traction-free surface, |$L_D^{*} \le h^{*}_F/e$| and is a function of dip angle, burial depth of the fault’s up-dip edge and friction coefficient. However, |$L_D^{*} \approx h^{*}_F/e$| for dip-slip faults where the up-dip edge of the fault is buried at a depth greater than or equal to |$h^{*}_F$|⁠. For sliding systems that can be treated as a thin layer, such as landslides, glaciers or ice streams, |$L_L^{*} = h^{*}_L/2$|⁠, where |$h^{*}_L$| is defined in Section 3.1. Conditions are also established for linear instability under velocity-strengthening friction on infinitely long faults that are parallel to a traction-free surface. Finally, since the focus of this paper is on linear stability, behaviour in the nonlinear regime (e.g. rupture localization or propagation) is not considered or examined; however simulations of the nonlinear governing equations are used to confirm some results of the linear analysis in Section 3.2.

2 METHODS

2.1 Linear stability analysis

Consider a fault of length L that obeys eq. (1) and either of eqs (2), and denote the position along the fault by |$\xi$|⁠. Assume also that the fault is embedded in a 2-D homogeneous elastic medium with shear modulus G and Poisson ratio |$\nu$|⁠. A linear stability analysis of the fault’s sliding motion can be conducted according to the following steps. (1) Write the system of nonlinear equations governing the evolution of sliding velocity |$v(\xi ,t)$| and state variable |$\theta (\xi ,t)$| along the fault. (2) Determine a uniform steady state of the system such that |$v(\xi ,t) = v_0$| and |$\theta (\xi ,t) = \theta _0$|⁠. Defining such a steady state requires that constant slip at rate |$v_0$| be imposed on extensions of the fault along the |$\xi$|-axis (i.e. backslip loading). (3) Obtain a linearized system of equations by computing the Jacobian matrix |$\boldsymbol {J}$| of the nonlinear system and evaluating it at the uniform steady state so that |$\boldsymbol {J}_0 = \boldsymbol {J}(v_0,\theta _0)$|⁠. (4) Determine the stability of the linear system by examining the eigenvalues of |$\boldsymbol {J}_0$|⁠. If any eigenvalue has a positive real part then the system is unstable.

Step 1. For quasi-static sliding, the velocity of the fault is governed by a balance between frictional resistance |$\tau _{F} = \mu \sigma$| and the shear stresses resolved upon the fault |$\tau = \tau _0 + \tau _E$|⁠, where |$\tau _E$| is the change in shear stress due to gradients in slip along the fault, and |$\tau _0$| is the shear stress on the fault in the absence of any slip. As with the shear stress, the normal stress on the fault is |$\sigma = \sigma _0 + \sigma _E$|⁠. The stress balance changes in time as |$\dot{\mu } \sigma = \dot{\tau }_E - \mu \dot{\sigma }_E$|⁠, and by making use of eq. (1), the sliding velocity of the fault can be written as

The evolution of the state variable can simply be written as |$\dot{\theta }(\xi , t) = H(v, \theta )$|⁠, since only the aging and slip laws are considered here and both state variable laws have the same linearization (e.g. Ruina 1983). The nonlinear governing equations for |$\dot{v}(\xi ,t)$| and |$\dot{\theta }(\xi ,t)$| are now represented by the functions |$F(v,\theta )$| and |$H(v, \theta )$|⁠.

When sliding is quasi-static, the changes in shear and normal stress are functions of the slip distribution |$\delta (\xi ,t)$|⁠, so that |${\tau }_E = T(\xi , \delta )$| and |$\sigma _E = N(\xi , \delta )$|⁠. The functions |$T(\xi , \delta )$| and |$N(\xi , \delta )$| are the stress change functions mentioned in the Introduction. These functions must be determined by solving the appropriate 2-D elasticity problem for a given fault geometry (see Appendix  C for example) and contain all necessary information about the elastic response of the system. These functions also have the property that |$\dot{\tau }_E = T(\xi , v)$| and |$\dot{\sigma }_E = N(\xi , v)$| (e.g. Viesca 2016a, b). Both |$T(\xi , \delta )$| and |$N(\xi , \delta )$| are equal to zero if there is no slip gradient.

Step 2. The uniform steady state of the system satisfies the conditions |$F(v_0, \theta _0) = 0$| and |$H(v_0, \theta _0) = 0$|⁠. Assume that the entire fault is sliding at steady state with velocity |$v_0$|⁠, such that |$T(\xi , v_0) = N(\xi , v_0) = 0$|⁠. This requires that the fault is subject to boundary conditions |$v(0,t) = v(L,t) = v_0$|⁠, and for nonzero |$v_0$| implies that constant sliding at |$v_0$| prevails on extensions of the fault along the |$\xi$|-axis (see dotted, red lines in Fig. 1). For both the aging and slip laws, |$H(v_0, \theta _0) = 0$| when |$\theta _0 = d_c/v_0$|⁠. These conditions satisfy |$F(v_0,\theta _0) = 0$|⁠, so the uniform steady state of the nonlinear system is |$(v_0, d_c/v_0)$|⁠.

Step 3. To linearize the equations about the uniform steady state, first define

where |$\boldsymbol {\mathrm{w}}(\xi ,t)$| is a small perturbation away from |$(v_0, \theta _0)$|⁠. Now the linearized equations can be written as |$\dot{\boldsymbol {\mathrm{w}}} = \boldsymbol {J}_0 \boldsymbol {\mathrm{w}}$|⁠. The Jacobian matrix |$\boldsymbol {J}_0$| is most conveniently expressed in terms of the dimensionless variables: |$\hat{t} = (v_0/d_c)t$|⁠, |$\hat{v} = v/v_o$| and |$\hat{\theta } = (v_0/d_c)\theta$|⁠, such that the uniform steady state becomes |$(\hat{v}_0, \hat{\theta }_0) = (1,1)$|⁠. Then |$\boldsymbol {J}_0$| can be written as

where |$\hat{F}_0 = \hat{F}(\hat{v}_0,\hat{\theta }_0)$|⁠, |$\hat{H}_0 = \hat{H}(\hat{v}_0,\hat{\theta }_0)$|⁠, |$\boldsymbol {I}$| is the identity matrix and |$\hat{T}_{\hat{v}}$| and |$\hat{N}_{\hat{v}}$| denote derivatives with respect to |$\hat{v}$|⁠. Some additional mathematical steps are provided in Appendix  A. The dimensions of |$\boldsymbol {J}_0$| will depend on whether the Jacobian is treated analytically or numerically.

Step 4. The eigenvalues and eigenvectors of |$\boldsymbol {J}_0$| determine solutions to the linearized system (⁠|$\dot{\boldsymbol {\mathrm{w}}} = \boldsymbol {J}_0 \boldsymbol {\mathrm{w}}$|⁠) of the form |$\boldsymbol {\mathrm{w}}(\xi ,t) \propto \boldsymbol {w}(k \xi ) e^{p t}$|⁠. The eigenvectors |$\boldsymbol {w}(k\xi )$| represent small spatial perturbations of wavenumber k to the uniform steady state. The eigenvalues p are the corresponding growth rates of those perturbations. If all of the eigenvalues of |$\boldsymbol {J}_0$| have a negative real part then the system is linearly stable. If any eigenvalue has a positive real part then the system is linearly unstable (e.g. Strogatz 2018).

2.1.1 Analytical stability analysis

For analytical results, |$\boldsymbol {I} \rightarrow 1$| in eq. (7) and |$\boldsymbol {J}_0$| can be treated as a |$2 \times 2$| matrix. Then the eigenvalues are found by solving the characteristic equation of |$\boldsymbol {J}_0$|⁠, such that

The eigenvalues do not need to be explicitly determined in cases where |$T_v$| and |$N_v$| are purely real functions. Instead, the stability of the system can be determined from conditions on |$\text{det}(\boldsymbol {J}_0) = -(b/a)\Gamma$| and |$\mathrm{Tr}(\boldsymbol {J}_0) = (b/a) (\Gamma + 1) - 1$|⁠, where |$\text{det}()$| and Tr() denote the determinant and trace of the matrix, respectively. The system is unstable if either |$\mathrm{Tr}(\boldsymbol {J}_0) > 0$|⁠, or |$\mathrm{det}(\boldsymbol {J}_0) < 0$| (e.g. Strogatz 2018, fig. 5.2.8); however, for all of the cases examined in this paper |$\mathrm{det}(\boldsymbol {J}_0) > 0$|⁠. The trace instability condition can be written as

Eq. (9) can be used to obtain analytical stability results. For example, in a spring-slider system |$N_v = 0$| and |$T_v = -K$|⁠, where K is the spring stiffness. Then solving eq. (9) for K will yield the usual relation for the critical spring stiffness (see Appendix B1).

Analytical results can be obtained for faults in a 2-D medium by specifying the functional form of the spatial perturbation |$\boldsymbol {w}(k \xi )$|⁠. For infinitely long faults there are no restrictions on the values of k because the fault has no boundaries. Then the general solution to the linear equations is

where |$A(k)$| is determined by a Fourier transform of the initial conditions |$\boldsymbol {\mathrm{w}}(\xi ,0)$| (Pivato 2010).

For a finite length fault with a uniform steady state velocity |$v_0$|⁠, the sliding velocity must remain |$v_0$| at the boundaries and so |$\boldsymbol {w}(k \xi )$| must be equal to zero at the boundaries. If the fault is defined over |$\xi = [0, L]$|⁠, the general solution to the linear equations that satisfies these boundary conditions is

where the constants |$\boldsymbol {A}_n$| are determined by Fourier series expansion of the initial conditions (Pivato 2010). The allowable wavenumbers are |$k = n\pi /L$| (for |$n = 1, 2, ...$|⁠), analogous to the normal modes on a vibrating string. Thus, the possible values of k are discrete on a finite length fault and scaled by the fault length L, and this has important consequences for the stability behaviour.

The general solutions for an initial perturbation given by eqs (10) and (11) also reveal why there is no condition on the spatial extent of a small perturbation for linearly unstable sliding. Since eqs (10) and (11) both represent a superposition of mode numbers, instability will occur for any initial perturbation |$\boldsymbol {\mathrm{w}}(\xi ,0)$| that contains wavelengths that have a positive growth rate p. For example, consider an infinitely long, velocity-weakening fault such that the critical perturbation wavelength is given by eq. (4). Any set of initial conditions |$\boldsymbol {\mathrm{w}}(\xi ,0)$| that contains wavelengths larger than |$h^{*}_F$| will have a positive growth rate and so be linearly unstable. Now consider a perturbation from steady sliding in the form of a finite delta function, such that |$\boldsymbol {\mathrm{w}}(\xi , 0) = \boldsymbol {\epsilon }$| at |$\xi = 0$|⁠, where |$\boldsymbol {\epsilon }$| is some small constant, and otherwise |$\boldsymbol {\mathrm{w}}(\xi ,0) = 0$|⁠. This set of initial conditions is unstable because it contains all wavelengths, since its Fourier transform is a constant. Since the finite delta function has no width, this means that there is no minimum slip patch length required to generate a linear instability. The same type of argument applies to finite faults as well and is further illustrated in Sections 3.1 and 3.2.

2.1.2 Numerical stability analysis

To conduct a linear stability analysis numerically, the fault can be discretized into |$n_e$| elements of length |${\rm d}\xi$|⁠. Then |$\boldsymbol {J}_0$| becomes a |$2n_e \times 2n_e$| block matrix where in general the upper-left block |${\partial {\hat{F}_0}} / {\partial {\hat{v}}}$| is dense and the other blocks are sparse diagonal matrices. In discrete form, the functions |$T(\xi ,v)$| and |$N(\xi ,v)$| become linear operators on the slip velocity. For example, |$T(\xi ,v) \rightarrow \sum _{j=1}^{n_e} T_{ij} v_j$| and |$T_v \rightarrow T_{ij}$|⁠, where |$T_{ij}$| is an |$n_e \times n_e$| matrix and |$v_j$| is a vector of length |$n_e$|⁠. All numerical results in this paper were obtained using a piecewise constant discretization of the stress change functions by assuming that slip is constant over regularly spaced elements along the fault.

The stability condition given by eq. (9) is only valid for |$2 \times 2$| matrices (e.g. Luís 2021). The eigenvalues must be explicitly calculated for numerical analysis (e.g. Viesca 2016a, b; Ray & Viesca 2017; Viesca 2023). The eigenvalues and eigenvectors can be directly computed using standard numerical routines; here the MATLAB functions eig and eigs are used. Numerically computing the eigenvalues will indicate if a fault system is stable or unstable for the specific set of RSF, elastic and geometrical parameters that define the system. Determining the conditions (if any) where stability changes requires an iterative assessment of the stability for different parameter values. In this paper, critical fault lengths are determined using a bisection method to locate values of L where the stability changes (to within |$\pm {\rm d}\xi /2$|⁠) while other properties are held constant.

2.2 Nonlinear simulations

In Section 3.2 the results of a limited set of simulations of the full nonlinear governing equations are presented to confirm some of the linear stability results. In these simulations the fault is loaded such that the steady state slip velocity along the entire fault is equal to |$v_0 = 10^{-9}$| m s⁻¹. These simulations use the aging law for state variable evolution and rather than eq. (1), use the regularized form of the rate and state friction equation (Rice & Ben-Zion 1996; Lapusta et al. 2000)

In discrete form the stressing rate balance at the centre of each fault element is

where |$\dot{\mu }$| is found from eq. (12) and |$\dot{\tau }_I$| is the radiation damping approximation for the inertial stressing rate (Rice 1993). The governing eq. (13) with eq. (12) and the aging law were solved along the entire length of the fault using a boundary element method implemented in MATLAB (see Data Availability statement for code availability).

3 RESULTS

Results are presented first for a finite length fault that is parallel to a traction-free surface, using a thin-layer approximation for the stress change functions. Analytical results are obtained for the thin-layer system, which provide insight into more complicated systems; results from numerical stability analysis for this system are presented as well. Next results are obtained numerically and semi-analytically for vertical strike-slip faults, as well as faults in an infinite full-space. Finally, dip-slip faults of any orientation in a system with a traction-free surface are examined.

3.1 Thin layer approximation

Consider a fault of length L that is parallel to a traction-free surface at a depth d (Fig. 1 B). In general, this system will have a non-zero |$N(\xi ,v)$| for in-plane sliding (see Section 3.3). However, when |$d \ll L_b = d_c G^{\prime }/(b \sigma _0)$| then |$N(\xi ,v)=0$| and the change in shear stress is (Viesca 2016a)

Note that eq. (14) is a special case of the stress change function for a dipping fault geometry illustrated in Fig. 1(C), as described in Section 3.3. The critical wavelength for an infinitely long fault in this system is (Viesca 2016b)

where |$L_{bh} = \sqrt{{\rm d} E^{\prime } d_c/(b\sigma _0)}$| (see Appendix B3.1).

Due to the simplicity of eq. (14), analytical results for the critical fault length and the wavelengths of unstable modes can be obtained for finite length faults in this system. By assuming a solution for |$v(\xi , t)$| of the form of eq. (11), the normalized derivative of the shear stress change function becomes |$\hat{T}_{\hat{v}}/b = -(n\pi L_{bh} / L)^2$| (see Appendix B3.2 for details). Then via eq. (9) the instability condition for the fault length becomes

Since the right hand side of eq. (16) is smallest at |$n=1$|⁠, the critical fault length is

Eq. (17) indicates that the fault becomes unstable when it is long enough that the wavelength |$\lambda$| of the first mode (⁠|$n=1$|⁠) of eq. (11) becomes equal to |$\lambda = 2\pi /k = 2L = h^{*}_L$|⁠.

The critical fault length |$L^{*}_L$| can also be numerically determined using the method described in Section 2.1.2. In this case by choosing a value of |$a/b$| then computing the stability of the system for different values of L. Then the critical fault length coincides with the value of L where the stability changes. Fig. 2(A) displays the results of this process for nine different values of |$a/b$| and shows that the numerically determined values of |$L^{*}_L$| agree with eq. (17).

As the fault length increases above |$L^{*}_L$|⁠, progressively higher mode numbers will become unstable and the wavelength of the highest unstable mode number will approach |$h^{*}_L$| as |$L \rightarrow \infty$|⁠. From eq. (16), the total number of unstable modes that a fault can host is |$n_T = \text{Fl}(2L/h^{*}_L)$|⁠, where |$\text{Fl}(q)$| gives the greatest integer less than or equal to some quantity q. The wavelength of the highest mode number |$n_T$| as a function of the fault length is

Eq. (18) predicts that |$\lambda _{n_T} \ge h^{*}_L$| and approaches |$h^{*}_L$| with a type of saw-tooth pattern as |$L \rightarrow \infty$| (Fig. 2B). This result can also be confirmed numerically by computing the wavelength of the eigenvector for the highest unstable mode as a function of the fault length for |$L > L^{*}_L$| (Fig. 2B). The close agreement between the analytical and numerical analyses both validates the numerical method and confirms the behaviour for finite length faults.

It is important to emphasize again that there is no minimum slip patch length required to generate an unstable sliding event for faults longer than |$L^{*}_L$|⁠. Since eq. (11) is a sum over all mode numbers, instability will occur for any set of initial conditions that contains a mode with a positive growth rate p. For example, consider a fault in the thin layer geometry with length |$L = L^{*}_L$|⁠, so that one unstable mode exists with wavelength |$\lambda = h^{*}_L$|⁠. Following the argument at the end of Section 2.1.1, an initial perturbation consisting of a finite delta function will be unstable, since the delta function can be represented as a sum over all mode numbers and therefore contains the unstable mode. This implies that there is no minimum perturbation length required to initiate a linear instability.

3.2 Vertical strike-slip faults in a half-space and full-space faults

Now consider a fault of length L embedded in a homogeneous full-space (Fig. 1A). For this system |$N(\xi ,v) = 0$| and the change in shear stress is given by (e.g. Segall 2010)

This stress change function is also valid for a vertical strike-slip fault in a half-space, in which case the integration is taken over |$[d, d+L]$| and |$G^{\prime } = G$| (Fig. 1C with |$\beta = 90^{\circ }$|⁠). Eq. (19) takes the form of a Hilbert transform for an infinitely long fault (⁠|$L \rightarrow \infty$|⁠). Then the critical wavelength |$h_F^{*}$| given by eq. (4) can be obtained from eq. (9) after applying a Fourier transform (see Appendix B2.1).

Analytical analysis using Fourier transforms cannot be applied to finite length faults due to the finite integration interval in eq. (19). However, the critical fault length can be obtained semi-analytically by using some results from Uenishi & Rice (2003) to solve eq. (9) (e.g. Ciardo & Viesca 2024). Details are provided in Appendix B2.2. The stability analysis can also be conducted numerically in the same manner as for the thin layer system, using the method described in Section 2.1.2. Both the semi-analytical and numerical stability analysis (Fig. 2A) show that the critical fault length for the full-space system is

Eq. (20) is an approximate equality in the absence of fully analytical results.

Eq. (20) is also supported by simulations of the full nonlinear governing equations following Section 2.2. The goal of these simulations is not to study nonlinear behaviour, but to confirm the linear stability results by using the full governing equations to simulate the growth or decay of an initial small perturbation from steady state. Fig. 3 shows the results of six sets of simulations using three values of |$a/b$| and two values of |$L_b$|⁠. Nine simulations, each with a different fault length, were run for each pair of (⁠|$a/b$|⁠, |$L_b$|⁠) values. In these simulations the initial conditions were set to the uniform steady state values, except for one element at the centre of the fault where |$v(\xi =0,t=0) = 0.99 v_0$|⁠. Hence the spatial extent of the initial perturbation is as small as the numerical discretization allows. Three additional sets of simulations for |$L_b = 1.6$| km were conducted with the perturbation applied to a single element at the edge of the fault.

The simulations were run until either consistent oscillations of maximum slip rate on the fault developed (i.e. a limit cycle), or the sliding velocity reached a uniform steady state such that |$v(\xi ,t) = v_0$|⁠. The critical fault length for each pair of (⁠|$a/b$|⁠, |$L_b$|⁠) values lies in the interval of fault lengths that separate growth and decay of the initial perturbation, as indicated by the maximum slip velocity. These critical fault lengths (normalized by |$h^{*}_F$|⁠) are |$0.37125 \pm 0.00125$| for |$a/b = 0.3, 0.7$| and |$0.37375 \pm 0.00125$| for |$a/b = 0.5$|⁠. These critical fault lengths are within |$2\, \rm per\,cent$| of the value given by eq. (20), and there is no dependence on the value of |$L_b$| or the location of the perturbation (Fig. 3).

Since the perturbation was restricted to a single fault element, the element length places an upper bound on any possible perturbation length scale required to initiate linear instability. The element length is set to |${\rm d}\xi = L_b/80$| in all of the simulations, and for the values of |$a/b$| used in Fig. 3 this corresponds to element lengths of |${\rm d}\xi \approx h^{*}_F/838$| to |$h^{*}_F/359$|⁠, since |$h^{*}_F = \pi L_b/(1 - a/b)$|⁠. These results further illustrate that |$h^{*}_F$| should not be interpreted as a minimum patch length needed to initiate a linearly unstable sliding event.

3.3 Dip-slip faults

Consider in-plane sliding on a fault that is dipping at an angle |$\beta$| relative to the traction-free surface of a homogeneous, elastic half-space (Fig. 1C). The up-dip edge of the fault is buried at a depth d below the traction-free surface. Both the full-space and parallel fault geometries are special cases of this dipping fault geometry. The full-space geometry is obtained when |$d \rightarrow \infty$|⁠, and the parallel fault geometry is obtained when |$d \ne 0$| and |$\beta = 0$|⁠.

Stress change functions for the half-space geometry are available in the literature (Dmowska & Kostrov 1973; Freund & Barnett 1976; Rudnicki & Wu 1995), and can be written as

where |$l = d / \sin {(\beta )}$|⁠, and |$\Psi (z,\beta )$| is an analytic function of the complex variable |$z = x + iy$| (England 2003). A similar expression holds for |$N(\xi ,v)$|⁠. A derivation of these functions is presented in Appendix  C. Note that these stress change functions are equivalent to using the Okada (1992) solutions for the middle of a very long dip-slip fault (e.g. Liu & Rice 2007).

3.3.1 Velocity-weakening behaviour

Critical fault lengths |$L^{*}_D$| for the dipping geometry can be determined by choosing a burial depth d and dip angle |$\beta$| and then conducting a numerical stability analysis as described in Section 2.1.2. Changing the value of L for fixed values of d and |$\beta$| corresponds to changing the down-dip depth of the fault. The stability calculation was repeated for dip angles in the range |$\beta = 0^\circ$| – |$90^\circ$| and burial depth values |$d/ h^{*}_F = 0$|⁠, |$10^{-3}$|⁠, |$10^{-2}$|⁠, |$10^{-1}$|⁠, 1 (the value |$d = 0$| was omitted for |$\beta = 0^\circ$|⁠). This process was carried out for values of |$\mu _0 = 0.2$|⁠, 0.6, 1, for both thrust and normal faults (Fig. 4).

The critical fault length |$L^{*}_D$| approaches the full-space value given by eq. (20) as |$d \rightarrow h^{*}_F$|⁠. Therefore |$L^{*}_D = L^{*}_F$| at depths |$d \ge h^{*}_F$| and Fig. 4 shows critical fault lengths for both thrust and normal faults in any possible orientation. For burial depths |$d < h^{*}_F$|⁠, the critical fault length is approximately a log-linear function of d (Fig. 5A).

The critical fault length |$L^{*}_D$| depends on the dip angle in a manner that is different for thrust and normal faults. There is also a secondary dependence on the value of |$\mu _0$| that depends on the sense of slip. For both thrust and normal faults, |$L^{*}_D$| increases with dip angle up to a value of |$20^\circ$| – |$40^\circ$|⁠, depending on the burial depth and sense of slip. For thrust faults, |$L^{*}_D$| then decreases to a secondary minimum before increasing again as |$\beta \rightarrow 90^\circ$|⁠. For normal faults, |$L^{*}_D$| reaches a maximum value then decreases as |$\beta \rightarrow 90^\circ$|⁠. Increasing the value of |$\mu _0$| decreases |$L^{*}_D$| for thrust faults, and does the opposite for normal faults. Values of |$L^{*}_D$| can become quite small on shallowly dipping faults that are near to the traction-free surface. In particular, as |$\beta \rightarrow 0^\circ$| on faults that break the surface (⁠|$d = 0$|⁠), |$L^{*}_D/h^{*}_F \rightarrow 10^{-2}$| on normal faults and appears to approach zero on thrust faults.

The dependence of |$L^{*}_D$| on |$\beta$| and |$\mu _0$| can mostly be explained by considering the on-fault stressing rates due to a uniform slip velocity distribution |$\Delta v$| on a dipping fault of length L with burial depth |$d = 0$|⁠. The elastic stressing rate on the fault is |$\dot{\tau }_E - \mu \dot{\sigma }_E$| (see Section 2.1, Step 1), which can be computed by evaluating the stress change functions at the centre of the fault |$\xi = L/2$| (e.g. Kato & Hirasawa 1997). The stressing rate has a dependence on |$\beta$| and |$\mu _0$| that shares some of the same features as that of |$L^{*}_D$|⁠; including similar behaviour as |$\beta \rightarrow 0^\circ$| and |$\beta \rightarrow 90^\circ$|⁠, and the same style of dependence on |$\mu _0$| for thrust and normal faults (Fig. 5B).

The stressing rate calculation also provides an explanation for why values of |$L^{*}_D$| become very small at shallow dip angles. Sliding instability develops when the frictional weakening rate |$\dot{\mu }$| is greater than the elastic stressing rate. The elastic stressing rate is approximately proportional to |$\beta /L$| for dip angles less than about |$10^\circ$| – |$20^\circ$| (Fig. 5B). Then for a given set of frictional parameters, when the dip angle is small only shorter length faults can relieve elastic stress faster than the frictional weakening rate. This leads to the results displayed in Fig. 4.

3.3.2 Velocity-strengthening behaviour

As noted in the Introduction, it is possible for unstable behaviour to occur on velocity-strengthening faults when a coupling between slip and normal stress exists, that is, when |$N(\xi ,v) \ne 0$|⁠. The parameter space for the dipping fault geometry is large; the stability behaviour can be expected to depend on frictional and elastic parameters |$\mu _0$|⁠, |$a/b$|⁠, |$L_b$|⁠; burial depth d; dip angle |$\beta$| and fault length L. Additionally, while normalization by |$h^{*}_F$| accounts for dependence on RSF and elastic parameters for velocity-weakening behaviour, |$h^{*}_F$| does not exist on velocity-strengthening faults. Therefore the results in this section are restricted to an infinitely long fault that is parallel to a traction-free surface, which reduces the parameter space to |$\mu _0$|⁠, |$a/b$| and a normalized burial depth |$d/L_b$|⁠. In this case eq. (8) can be used to determine the stability of the system (see Appendix B4 for stress change functions and details).

Fig. 6 shows the results of choosing values of |$\mu _0$| and |$d/L_b$|⁠, then determining the maximum value of |$(a/b)_c = (a/b) > 1$| that satisfies |$\text{Re}(p) > 0$| in eq. (8). One striking feature of the results is that unstable behaviour only exists at depths greater than some minimum value that is very well approximated by

The details of obtaining eq. (22) are provided in Appendix B4. At shallower depths there are no unstable solutions to eq. (8) for |$(a/b) > 1$|⁠. This shallow, stable region is not related to the thin layer limit that occurs at |$d/L_b \ll 1$|⁠. Where unstable behaviour occurs, for constant |$\mu _0$| there is depth at which |$(a/b)_c$| reaches a maximum value. While for constant |$d/L_b$|⁠, values of |$(a/b)_c$| increase monotonically with |$\mu _0$|⁠, so the velocity-strengthening instability is enhanced when friction is higher.

An extensive parameter study to determine the effects of finite fault length and dip angle is beyond the scope of this study. However, some insight can be gained by examining the critical wavelengths that correspond to the values of |$(a/b)_c$|⁠. Each value of |$(a/b)_c$| shown in Fig. 6(A) occurs at some critical wavenumber |$k_c$| that is in the neighbourhood of |$k_c/d \approx 1$| regardless of the value of |$(a/b)_c$| (Fig. 6B). By analogy with the velocity-weakening results, if |$L^{*} \approx h^{*}/e = 2\pi / (e k_c)$|⁠, then for any set of values |$[\mu _0, (a/b)_c, d/L_b]$| taken from Fig. 6(A), the fault length would have to be |$L \ge 2 \pi d /e$| for unstable behaviour to occur.

4 DISCUSSION

4.1 Some theoretical considerations

A main result in this paper is that the linear stability of frictional sliding depends on overall fault length. The critical fault length |$L^{*}$| for finite length faults replaces the concept of the critical wavelength |$h^{*}$| for infinitely long faults. For vertical strike-slip faults or faults in a full-space |$L^{*} \approx h^{*}_F/e$| (Section 3.2), while for dip-slip faults the critical fault length is a function of the dip angle and burial depth (Section 3.3). For sliding systems that can be treated as a thin layer, |$L^{*} = h^{*}_L/2$| (e.g. landslides, glaciers or ice streams; Section 3.1). Velocity-weakening faults are linearly unstable if they are longer than |$L^{*}$|⁠, in the sense that an infinitesimally small perturbation to the steady state will grow exponentially until the linear approximation is no longer valid. The results also illustrate that there is no minimum length scale associated with initial perturbations, leading to the conclusion that the critical patch length |$l_c$| obtained from spring-slider models does not have any relevance to the linear stability behaviour of finite or infinite length faults in an elastic continuum. This may explain why studies of length scales related to nonlinear stability and nucleation behaviour have not found any connection with |$l_c$| (e.g. Rubin & Ampuero 2005; Ampuero & Rubin 2008).

The concept of a critical fault length also leads to a new definition of conditional stability. Velocity-weakening faults that are shorter than |$L^{*}$| are conditionally stable in that they are stable in the linear regime, but large perturbations out of the linear regime could generate unstable sliding events (e.g. Gu et al. 1984). In the classical model of fault stability, regions of a fault are assigned stability properties (stable, unstable or conditionally stable) based on the linear stability analysis of a spring-slider system and variations in |$(a-b)$| along the fault (Scholz 1998, fig. 2). The results in this paper make use of fault systems with constant values of |$(a-b)$|⁠; however the numerical method for linear stability analysis can easily be applied to systems where the frictional properties vary along the fault. Thus, it is possible to obtain linear stability results for heterogeneous fault systems as a whole, rather than apply results from spring-slider systems to individual sections of a fault. The linear stability of such fault systems can be expected to depend on geometrical aspects (e.g. dip angle, burial depth) as well as the frictional properties in both the velocity-weakening and velocity-strengthening portions of the fault (e.g. Skarbek et al. 2012; Dublanchet et al. 2013; Ray & Viesca 2017; Yabe & Ide 2017; Luo & Ampuero 2018).

4.2 Influence of the free surface and burial depth

Proximity to a traction-free surface, as measured by |$h^{*}_F$| or |$L_b$|⁠, has a significant influence on stability properties. Since both |$h^{*}_F$| and |$L_b$| are inversely proportional to effective normal stress, the normalized burial depths in Figs 4 and 6 are smaller on faults with high pore fluid pressure. This means that the influence of the free surface is enhanced on overpressured fault systems. High pore pressure leads to smaller normalized critical fault lengths, but larger values of |$h^{*}_F$|⁠. If the burial depth is less than |$h^{*}_F$|⁠, then the free surface will influence the stability behaviour. This effect should be important in the shallow regions of subduction zones and in areas of induced seismicity where pore pressures can be elevated. Particularly on subduction megathrust plate boundaries, the combination of shallow dip angles and high pore pressures should lead to very small normalized critical fault lengths.

The effect of shallow burial depth on unstable behaviour for velocity-strengthening faults is more complicated. A fault parallel to the free surface should be the most unstable geometry for a non-zero burial depth d, since on a dipping fault the depth from the traction-free surface will increase with down-dip distance. The values of |$(a/b)_c$| for the infinite fault system in Fig. 6(A) are close to velocity-neutral, so it seems reasonable to assume that values of |$(a/b)_c$| would be even closer to unity on finite length, dipping faults that are buried. However, the velocity-weakening results show that intersecting the free surface causes a significant reduction in stability; |$L^{*}_D/h^{*}_F$| decreases logarithmically with decreasing |$d/h^{*}_F$|⁠. So it is possible that values of |$(a/b)_c$| may be larger on dipping faults where |$d = 0$|⁠. Certainly more work is needed to understand this behaviour.

Multiple effects have been described that can cause unstable sliding on velocity-strengthening faults: contrasting elastic parameters across a fault (Rice et al. 2001; Ranjith 2014); poroelasticity (Heimisson et al. 2019); ‘fault valve’ behaviour (Ozawa et al. 2024); and proximity to a traction-free surface (this paper; Aldam et al. 2016). All of these features are commonplace in fault systems as well as in other frictional systems like landslides and ice streams. For example, all of these effects could be present in the shallow regions of subduction zones and may contribute towards enabling shallow slow-slip events (e.g. Saffer & Wallace 2015), or influencing the behaviour of tsunami earthquakes (e.g. Bilek & Lay 2002).

5 CONCLUSION

The results in this paper show how even simple types of geometrical complexity can change stability behaviour. Using numerical methods makes it possible to conduct linear stability analyses for a wide range of fault systems that cannot be examined using analytical techniques. Some examples of systems for which stress change functions are available in the literature are multifault systems and non-planer faults in a 3-D homogeneous elastic half-space (Okada 1992; Meade 2007). Functions are also available for different types of viscoelastic geometries (e.g. Segall 2010; Lambert & Barbot 2016). Heterogeneous on-fault frictional properties can be used with any existing stress change functions (e.g. Ray & Viesca 2017). Finally, numerical stability methods could also be extended to include dilatancy and changes in pore pressure, or other types of frictional constitutive behaviour (e.g. Segall & Rice 1995; Chen & Spiers 2016; Barbot 2022).

ACKNOWLEDGEMENTS

This work was supported by the National Science Foundation under Grant No. EAR-2245540. The author thanks D. Saffer, H. Savage and R. Viesca for discussions that helped and influenced this paper. Colours for Fig. 6 were generated using a perceptually uniform colour map (Crameri 2019; Crameri et al. 2020). The author also thanks Sohom Ray, the editor Shiqing Xu and an anonymous reviewer for comments that improved the manuscript.

DATA AVAILABILITY

All of the calculations and figures in this paper can be reproduced using a MATLAB package RSFaultZ available at https://github.com/rmskarbek/RSFaultZ (Skarbek 2024). The m-files for automatically generating figures are stored in the github repository directory: RSFaultZ/examples/stability.

REFERENCES

APPENDIX A: LINEARIZATION OF RSF EQUATIONS

Additional mathematical details are provided here for obtaining the Jacobian matrix given by eq. (7). First, using eqs (2) and (5), the linearized equations can be written as

and

Eqs (A1) and (A2) can be used to define a dimensional Jacobian. The elements of eq. (7) are obtained after changing to the dimensionless variables defined by |$\hat{t} = (v_0/d_c)t$|⁠, |$\hat{v} = v/v_o$| and |$\hat{\theta } = (v_0/d_c)\theta$|⁠. Dimensionless stress change functions are obtained by normalizing stresses by |$\sigma _0$|⁠. So for example, |$\dot{\tau }_E = T(\xi ,v) = (\sigma _0 v_0 / d_c) \hat{T}$|⁠.

APPENDIX B: ANALYTICAL LINEAR STABILITY RESULTS

B1 Spring-slider

The shear stress change in the basic spring-slider model is

where K is a normalized spring stiffness with units of [Stress / Length]. Using the same dimensionless variables defined in  A, the dimensionless shear stress change function is

Inserting the derivative of eq. (B2) with respect to |$\hat{v}$| into eq. (9) and setting the left-hand side equal to zero yields the critical stiffness |$K_c = \sigma _0(b-a)/d_c$|⁠.

B2 Full-space

B2.1 Infinite fault

For infinite faults the critical wavelength can be found by searching for solutions of the form |$v(\xi ,t) = A \exp (pt + i k \xi )$|⁠. For a full-space, the shear stress change function can be found by substituting this expression into eq. (19), for |$L \rightarrow \infty$|⁠; this is essentially the method used by (Rice et al. 2001):

After making a change of variables |$u = s - \xi$|⁠, eq. (B3) becomes

where the integral in the first line is a Fourier Transform of |$1/u$| and is equal to |$i \pi \text{sgn}(k)$|⁠. Using the previously defined dimensionless variables, but leaving k in dimensional form, the critical wavenumber |$k_c$| from eq. (9) is

which leads to eq. (4) since the critical wavelength is defined as |$h^{*}_F = \lambda _c = 2\pi / k_c$|⁠.

B2.2 Finite fault

The critical fault length for a finite fault in a full-space can be found by first normalizing eq. (19) and inserting the result into eq. (9). Again using the previously defined dimensionless variables and also normalizing lengths by |$L/2$|⁠, eq. (19) becomes

Substituting this result into eq. (9) and rearranging terms yields

Finally, integrating eq. (B7) with respect to |$\hat{v}$| provides an equation for the critical fault length

Note that the negative sign on the right-hand side of eq. (B7) has been absorbed into the denominator of the integral argument in eq. (B8).

Eq. (B8) is in the form of an eigenvalue equation where |$\pi L / h^{*}_F$| is the eigenvalue and |$\hat{v}(\hat{\xi })$| is the corresponding eigenvector. Note that this eigenvalue equation is distinct from what is described in Section 2.1. The critical fault length is found from the smallest eigenvalue of eq. (B8)

The numerical factor in eq. (B9) has previously been determined by Uenishi & Rice (2003); see eqs (12) and (13) and Appendix B1 in that paper (see also Appendix B in Ciardo & Viesca 2024). Finally, solving eq. (B9) for L provides the critical fault length given by eq. (20).

B3 Thin layer

B3.1 Infinite fault

The critical wavelength for the thin layer system can be found by following the same procedure for the full-space system, but using eq. (14) for the shear stress change function.

Using the dimensionless variables as before, eq. (9) becomes

with |$L_{bh} = \sqrt{{\rm d} E^{\prime }d_c/(b\sigma _0)}$| (Viesca 2016b). Solving eq. (B11) for the critical wavelength leads to eq. (15) for |$h^{*}_L$|⁠.

B3.2 Finite fault

For a finite fault the deviation of the sliding velocity from steady state takes the form |$v(\xi ,t) - v_0 = a_n e^{p t} \sin (n \pi \xi /L)$|⁠. After substituting this into eq. (14), the shear stress change becomes

such that

Eq. (B13) can be normalized using the dimensionless quantities defined in Appendix  A and remembering that T has units of [stress/time], then

as in Section 3.1. Finally, the critical fault length is obtained by substituting eq. (B14) into eq. (9), which yields eq. (16).

B4 Velocity-strengthening layer

The stress change functions for in-plane sliding on an infinitely long fault that is parallel to a traction-free surface at a depth d are (e.g. Viesca 2016a)

and

The stability of this system is most easily determined after applying a Fourier transform. Using the Fourier transform pair:

eqs (B15) and (B16) become

and

where tildes denote transformed quantities. Note that these functions are provided by Viesca (2016a) using a different transform pair.

The eigenvalues p can then be computed from eq. (8) after defining |$\Gamma = (1/b)(\widetilde{T}_{\tilde{v}} - \mu _0 \widetilde{N}_{\tilde{v}})$| using eqs (B19) and (B20). Using the dimensionless variables, and also defining |$\hat{k} = d k$| yields

The resulting equation for p is complex and depends on the values of |$(a/b)$|⁠, |$(L_b/d)$|⁠, |$\mu _0$| and the dimensionless wavenumber |$\hat{k}$|⁠. The results in Fig. 6 were obtained through an iterative process by solving for p numerically as a function of |$\hat{k}$| for chosen values of |$L_b/d$| and |$\mu _0$|⁠. For each pair of values |$(L_b/d, \mu _0)$|⁠, |$p(\hat{k})$| was first determined for a value of |$(a/b) < 1$|⁠, which guarantees that |$\text{Re}[p(\hat{k})] > 0$| for some value of |$\hat{k}$|⁠; numerical tests showed that the maximum value of |$p(\hat{k})$| occurs in the vicinity of |$\hat{k} \approx 1$|⁠. This process was then repeated for incrementally larger values of |$(a/b)$| until |$\text{Re}[p(\hat{k})] < 0$| for all values of |$\hat{k}$|⁠, which determines the values of |$(a/b)_c$| shown in Fig. 6(A).

The minimum depth for unstable behaviour can be approximately determined by solving for p for a specific value of |$\hat{k}$|⁠. From the results in Fig. 6(B), the stability boundary occurs at |$\hat{k} \approx 0.5$|⁠, so that

Additionally, that stability boundary occurs at |$(a/b) = 1$|⁠. With these values of |$\hat{k}$| and |$(a/b)$|⁠, eq. (8) becomes

Now p can be solved for using a procedure described in Rice et al. (2001). First, Fig. 6 indicates that for a constant value of |$\mu _0$|⁠, the real part of p changes sign as |$d/L_b$| increases from zero. The sign change occurs at |$p = i \rho$|⁠; substituting this into eq. (B23) yields

Eq. (B24) is satisfied when both its real and imaginary parts are equal to zero. Setting the real part equal to zero provides an equation for |$\rho$| in terms of |$\mu _0$| and |$(L_b/d)$|⁠:

Finally, inserting eq. (B25) into the imaginary part of eq. (B24) and setting it equal to zero provides an equation for |$d/L_b$| as a function of |$\mu _0$|⁠. The best way to execute this final step is using a symbolic math program. The solution is

APPENDIX C: DIP-SLIP FAULTS

Consider an edge dislocation in a 2-D homogeneous elastic body. The dislocation induces displacement and stress fields throughout the elastic body that can be represented in terms of two complex potentials, |$\omega (z)$| and |$\Omega (z)$|⁠, that are analytic functions of z (e.g. England 2003; Bower 2009). The complex coordinate z is defined as |$z = x + i y = r e^{i \phi }$| where |$(r,\phi )$| are radial coordinates with |$\phi$| measured from the x-axis in the direction of the y-axis.

For the dipping fault system shown in Fig. 1(C), the fault is located at |$\beta _0 = \pi - \beta$| along |$l \le r \le l + L$|⁠, where |$l = d / \sin (\beta )$| (also note that |$\xi = r$|⁠). The stress change functions can be obtained by considering a distribution of dislocations along the fault, and computing the shear and normal stresses that these dislocations induce on the fault itself. The first and most important step is to determine the complex potentials for a single dislocation placed at |$z_0 = r_0 e^{i \beta _0}$|⁠, with Burger’s vector |$b e^{i \beta _0} = b \cos (\beta _0) + i b \sin (\beta _0)$| (e.g. Freund & Barnett 1976).

In the x-y plane the stress and displacement fields are given by:

where primes denote derivatives with respect to z, and bars denote complex conjugates (e.g. section 2.5 in England 2003). The displacements are denoted by |$u_x$|⁠, |$u_y$|⁠; the normal stresses by |$\sigma _x$| and |$\sigma _y$| and |$\sigma _{xy}$| is the shear stress. The normal and shear stresses on the fault can be obtained in the radial coordinate system, in which case the stresses are

For a half-space with a traction-free surface at |$y = 0$|⁠, |$z = x$|⁠, the potentials can be written as

where |$\Omega _0(z)$| and |$\omega _0(z)$| are the potentials for a full-space, and so will produce tractions along |$z = x$|⁠; while |$\Omega _1(z)$| and |$\omega _1(z)$| are additional potentials that clear the tractions along |$z = x$|⁠. The full-space potentials are given by (e.g. Bower 2009, section 5.3.12)

where

The additional potentials can be found using a variety of methods (e.g. Dmowska & Kostrov 1973; Freund & Barnett 1976). Here, the additional potentials are computed using the process of analytic continuation (e.g. section 3.5 in England 2003), and are given by

Substituting these definitions for |$\Omega _1(z)$| and |$\omega _1(z)$| into eqs (C7) along with the results for |$\Omega _0(z)$| and |$\omega _0(z)$|⁠, the potentials for an edge dislocation in a half-space are:

Note that eqs (A4) and (A5) in Rudnicki & Wu (1995) are the derivatives of eqs (C13) and (C14).

The normal |$\sigma _\phi$| and shear |$\sigma _{r \phi }$| stresses on the fault due to a single dislocation are given by the real and imaginary parts of eq. (C5), evaluated using eqs (C13) and (C14) at values of z corresponding to |$\phi = \beta _0$| and |$l \le \xi \le l + L$|⁠. For a distribution of dislocations along the length of the fault, the resultant Burger’s vector between neighbouring points |$\xi$| and |$\xi + {\rm d}\xi$| is |$b = (\partial {\delta }/ \partial {\xi }) {\rm d}\xi$|⁠, where |$\delta (\xi )$| is slip on the fault (Freund & Barnett 1976; Weertman 1996). The stress change functions are found by integrating over the length of the fault, such that

where the potentials are evaluated using eqs (C13) and (C14) at |$z = \xi e^{i \beta _0}$| and |$z_0 = s e^{i \beta _0}$|⁠. Finally, note that it is possible to write the integrands in eqs (C15) and (C16) explicitly in terms of |$\xi$| and |$\beta _0$|⁠, however the resulting expressions are extremely cumbersome [see for example eqs (13) in Freund & Barnett (1976); eqs (3.1)–(3.2) in Dmowska & Kostrov (1973); or eqs (A6)–(A11L) in Rudnicki & Wu (1995)]. For numerical computations it is most concise to compute the stresses using the individual equations listed above.