A computational method for earthquake cycles within anisotropic media

SUMMARY

We present a numerical method for the simulation of earthquake cycles on a 1-D fault interface embedded in a 2-D homogeneous, anisotropic elastic solid. The fault is governed by an experimentally motivated friction law known as rate-and-state friction which furnishes a set of ordinary differential equations which couple the interface to the surrounding volume. Time enters the problem through the evolution of the ordinary differential equations along the fault and provides boundary conditions for the volume, which is governed by quasi-static elasticity. We develop a time-stepping method which accounts for the interface/volume coupling and requires solving an elliptic partial differential equation for the volume response at each time step. The 2-D volume is discretized with a second-order accurate finite difference method satisfying the summation-by-parts property, with boundary and fault interface conditions enforced weakly. This framework leads to a provably stable semi-discretization. To mimic slow tectonic loading, the remote side-boundaries are displaced at a slow rate, which eventually leads to earthquake nucleation at the fault. Time stepping is based on an adaptive, fourth-order Runge–Kutta method and captures the highly varying timescales present. The method is verified with convergence tests for both the orthotropic and fully anisotropic cases. An initial parameter study reveals regions of parameter space where the systems experience a bifurcation from period one to period two behaviour. Additionally, we find that anisotropy influences the recurrence interval between earthquakes, as well as the emergence of aseismic transients and the nucleation zone size and depth of earthquakes.

Creep and deformation, Numerical approximations and analysis, Numerical modelling, Seismic anisotropy

1 INTRODUCTION

Modelling the full earthquake cycle poses numerous computational challenges. Interseismic periods between fault rupture last hundreds of years, punctuated by earthquakes that evolve on a timescale of seconds. The spatial scales that must be considered also encompass many orders of magnitude. Fault length is measured in kilometres while the process zone, an area directly behind the tip of a propagating rupture, must be resolved on the order of millimetres when using laboratory-measured parameters (Noda et al. 2009). Additionally, faults in nature have nonplanar geometries, and the physical makeup of the materials that surround earthquake faults is complex and varying. Material anisotropy is present in the Earth’s crust, the upper mantle, the transition zone, the D” layer and the inner core (Long & Becker 2010), and seismic anisotropy can be observed through shear wave splitting, that is, when a shear wave splits into two components with different propagation speeds. This splitting has been observed in most igneous, metamorphic and sedimentary rocks in the Earth’s crust (Crampin & Lovell 1991) and is used to measure anisotropy along fault rupture zones (Cochran et al. 2003). While seismic anisotropy is present in the real world, many cycle models make the simplifying assumption of isotropic material properties (Lapusta et al. 2000; Erickson & Dunham 2014; Allison & Dunham 2018).

Methods currently used for simulating earthquake cycles can generally be broken into two broad categories: spectral boundary integral techniques, and numerical discretizations of off-fault volumes like finite difference and finite element methods. Hajarolasvadi & Elbanna (2017) detail the benefits and drawbacks of spectral boundary integral methods, namely, that they are computationally efficient [as they reduce 2-D problems to 1-D and 3-D problems to 2-D (Cochard & Madariaga 1994; Geubelle & Rice 1995)] and require no artificial truncation of the computational domain. However, these methods are currently limited to the simplifying assumption that the Earth’s material properties are homogeneous, isotropic and linear elastic. This motivates the use of volume-based numerical methods, such as finite element and finite difference methods which can account for material anisotropy, heterogeneity and off-fault plasticity (Kaneko et al. 2008; Erickson & Dunham 2014; Erickson et al. 2017; Allison & Dunham 2018).

In this work, where our focus is on anisotropic material properties, we elect to use a finite difference formulation satisfying a summation-by-parts (SBP) rule, with weak enforcement of boundary conditions through the simultaneous-approximation-term (SAT), which have the desirable property that the discretization is provably energy stable. This computational framework is an extension of that of Erickson & Dunham (2014) to incorporate material anisotropy. Our main focus is a parameter exploration of the homogeneous, anisotropic problem. The paper is organized as follows: In Section 2 we define the governing equation and constitutive relations for an anisotropic elastic material. In Section 3 we provide details of the spatial discretization and derive conditions that render the semi-discrete equations stable. In Section 4 we provide details of the frictional fault that forms one boundary of the domain and describe the adaptive Runge–Kutta-based time-stepping method. Section 5 is a linear stability analysis of frictional sliding for the anisotropic problem that extends the analysis done in Ranjith & Gao (2007). To verify our computational strategy, we perform convergence tests in Section 6 and confirm that our numerical solution is converging at the expected rate. Results from our parameter varying study are detailed in Section 7.

2 GOVERNING EQUATIONS

We consider a vertical, strike-slip fault embedded in a 2-D volume given by (y, z) ∈ ( − L_y, L_y) × (0, L_z). We assume that the only non-zero component of displacement, denoted as u, occurs in the x-direction and that motion is invariant in this direction, so that u = u(t, y, z). The fault, at y = 0, serves as a frictional interface embedded in an anisotropic, homogeneous material. Across the fault interface the components of traction are taken to be equal in magnitude but opposite in sign and the displacement u is allowed to have a jump. The jumps in displacement are governed by a friction law which couples the volume solution to a set of local auxiliary state variables governed by an ODE; the details of the frictional framework are given in Section 4.

Since the domain is symmetric about the fault and the material properties are homogeneous, the solution across the fault will be antisymmetric. Exploiting this symmetry allows us to consider the one sided domain (y, z) ∈ Ω = (0, L_y) × (0, L_z). In this setting, we extend the method from Erickson & Dunham (2014) to the anisotropic case by substituting Hooke’s law into momentum balance. The resulting 2-D anisotropic wave equation,

$$ \begin{eqnarray*} \rho \displaystyle\frac{\partial ^{2} u}{\partial t^{2}} &=& \displaystyle\frac{\partial }{\partial y}\left[\mu _{1}\displaystyle\frac{\partial u}{\partial y} + \mu _2 \displaystyle\frac{\partial u}{ \partial z}\right] \nonumber \\ &&+\displaystyle\frac{\partial }{\partial z}\left[\mu _{2}\displaystyle\frac{\partial u}{\partial y} + \mu _3 \displaystyle\frac{\partial u}{ \partial z}\right],\quad (y,z)\in \Omega , \end{eqnarray*}$$

(1)

contains mixed derivative terms, unlike its isotropic counterpart, due to the more complex form of the anisotropic material stiffness tensor. Here, ρ is the material density and μ₁, μ₂ and μ₃ are the elastic moduli; we assume that these satisfy μ₁ > 0, μ₃ > 0 and |$\mu _1 \mu _3 \gt \mu _2^2$|⁠.

The relevant components of traction that will be needed later are

$$ \begin{eqnarray*} \sigma _{xy}(t,y,z) &= \mu _1 \displaystyle\frac{\partial u}{\partial y} + {\mu} _2 \displaystyle\frac{\partial u}{\partial z} \end{eqnarray*}$$

(2a)

$$ \begin{eqnarray*} \sigma _{xz}(t,y,z) &= \mu _2 \displaystyle\frac{\partial u}{\partial y} + \mu _3 \displaystyle\frac{\partial u}{\partial z}. \end{eqnarray*}$$

(2b)

We impose the following boundary conditions on ∂Ω, as in Erickson & Dunham (2014), noting that in the anisotropic setting σ_xy and σ_xz depend on both y and z, unlike in the isotropic case, where σ_xy only depends on y and σ_xz only depends on z:

$$ \begin{eqnarray*} u(t,0,z) &= g_L(t,z), \end{eqnarray*}$$

(3a)

$$ \begin{eqnarray*} u(t,L_y,z) &= g_R(t,z), \end{eqnarray*}$$

(3b)

$$ \begin{eqnarray*} \sigma _{xz}(t,y,0) &= 0, \end{eqnarray*}$$

(3c)

$$ \begin{eqnarray*} \sigma _{xz}(t,y,L_z) &= 0. \end{eqnarray*}$$

(3d)

Condition (3c) corresponds to the Earth’s free surface and condition (3d) corresponds to the assumption that the material below depth L_z exerts no traction on the overlying material. The displacement boundary condition g_L(t, z) is determined by the friction law and g_R(t, z) imposes the remote tectonic loading; see Section 4. In order to derive parameters in the discretization, in eq. (1) we have retained the inertial term ∂²u/∂t². Later, in order to make the problem more computationally tractable, this inertial term will be replaced with the radiation damping approximation (Rice 1993) which will result in a modification to boundary condition (3a).

2.1 Energy-boundedness of the solution

Although Erickson & Dunham (2014) was a provably stable formulation, the additional mixed derivative terms in eq. (1) and additional terms present in eqs (2a) and (2b) necessitate a new analysis. As such, to ensure that the initial boundary value problem (1)–(3) is well-posed we use the energy method. We assume homogeneous boundary conditions, with the understanding that the analysis for zero boundary data can be extended to non-homogeneous boundary conditions via Duhamel’s principle. Letting || · || denote the L² norm, multiplying eq. (1) by |$\dfrac{\partial u}{\partial t}$|⁠, integrating over Ω and applying Green’s theorem on the right-hand side yield

$$ \begin{eqnarray*} \dfrac{1}{2}\displaystyle\frac{\partial }{\partial t} \left\Vert \rho \displaystyle\frac{\partial u}{\partial t} \right\Vert ^{2} &=& B_r +B_s +B_d + B_f - \dfrac{1}{2}\displaystyle\frac{\partial }{\partial t} \iint \limits _{\scriptstyle 0\, 0}^{\scriptscriptstyle L_y L_z} {\begin{bmatrix}\dfrac{\partial u}{\partial y} & \dfrac{\partial u}{\partial z} \end{bmatrix}}\nonumber \\ && \times {\begin{bmatrix}\mu _1 & \mu _2\\ \mu _2 & \mu _3 \end{bmatrix}} {\begin{bmatrix}\dfrac{\partial u}{\partial y} \\ \dfrac{\partial u}{\partial z} \\ \end{bmatrix}} \, \mathrm{ d}z\, \mathrm{ d}y, \end{eqnarray*}$$

(4)

where the boundary terms are given by

$$ \begin{eqnarray*} B_r &=& \int \limits _{0}^{L_{z}} \left( \mu _1 \dfrac{\partial u}{\partial t} \dfrac{\partial u}{\partial y} + \mu _2 \dfrac{\partial u}{\partial t} \dfrac{\partial u}{\partial y} \right) \Big \vert _{y=L_y} \nonumber \\ \mathrm{ d}z &=& \int \limits _{0}^{L_{z}} \left(\dfrac{\partial u}{\partial t}\sigma _{xy}(t,y,z) \right)\Big \vert _{y=L_y} \, \, \mathrm{ d}z, \end{eqnarray*}$$

(5a)

$$ \begin{eqnarray*} B_f &=& -\int \limits _{0}^{L_{z}} \left(\mu _1 \dfrac{\partial u}{\partial t} \dfrac{\partial u}{\partial y} + \mu _2 \dfrac{\partial u}{\partial t} \dfrac{\partial u}{\partial y} \right) \Big \vert _{y=0} \nonumber \\ \mathrm{ d}z &=& \int \limits _{0}^{L_{z}} \left(\dfrac{\partial u}{\partial t}\sigma _{xy}(t,y,z) \right)\Big \vert _{y=0} \, \, \mathrm{ d}z, \end{eqnarray*}$$

(5b)

$$ \begin{eqnarray*} B_d &=& \int \limits _{0}^{L_{y}} \left(\mu _2 \dfrac{\partial u}{\partial t} \dfrac{\partial u}{\partial z} + \int \limits _{0}^{L_{y}} \mu _3 \dfrac{\partial u}{\partial t}\dfrac{\partial u}{\partial z} \right) \Big \vert _{z=L_z}\nonumber \\ \mathrm{ d}y &=& \int \limits _{0}^{L_{z}}\left( \dfrac{\partial u}{\partial t}\sigma _{xz}(t,y,z) \right)\Big \vert _{z=L_z} \, \, \mathrm{ d}y, \end{eqnarray*}$$

(5c)

$$ \begin{eqnarray*} B_s &=& -\int \limits _{0}^{L_{y}} \left( \mu _2 \dfrac{\partial u}{\partial t} \dfrac{\partial u}{\partial z}+ \int \limits _{0}^{L_{y}} \mu _3 \dfrac{\partial u}{\partial t}\dfrac{\partial u}{\partial z} \right) \Big \vert _{z=0}\nonumber \\ \mathrm{ d}y &=& \int \limits _{0}^{L_{z}} \left(\dfrac{\partial u}{\partial t}\sigma _{xy}(t,y,z) \right)\Big \vert _{z=0} \, \, \mathrm{ d}y. \end{eqnarray*}$$

(5d)

Letting |$\mathbf {U} = [u_1 \quad u_2]^T$| allows us to define the norm

$$ \begin{eqnarray*} \left\Vert \mathbf {U}\right\Vert ^2_{M_{\mu }}=\iint \limits _{\scriptstyle 0\, 0}^{\scriptscriptstyle L_z L_y}{ {\mathbf {U}^{T}\mathbf {M}_{\mu }\mathbf {U}}} \, \mathrm{ d}y\, \mathrm{ d}z, \quad \mathbf {M}_{\mu }= {\begin{bmatrix}\mu _1 & \mu _2\\ \mu _2 & \mu _3 \end{bmatrix}}. \end{eqnarray*}$$

(6)

Here, |$\mathbf {M}_{\mu }$| is symmetric positive-definite due to the restrictions on the shear moduli given after eq. (1). The energy method is now complete, since we can write (4) as

$$ \begin{eqnarray*} \dfrac{1}{2}\dfrac{\partial }{\partial t} \left(\left\Vert \sqrt{\rho } \dfrac{\partial u}{\partial t}\right\Vert ^{2} + \left\Vert {\begin{bmatrix}\dfrac{\partial u}{\partial y} \dfrac{\partial u}{\partial z} \end{bmatrix}}^T\right\Vert _{M_{\mu }}^{2}\right)= B_r +B_s +B_d + B_f, \nonumber \\ \end{eqnarray*}$$

(7)

where the terms on the left-hand side of eq. (7) correspond to the rate of change of kinetic and strain energy, respectively, and the terms on the right-hand side represent rate of work done on the elastic body by tractions on the boundaries. For zero boundary data, B_r, B_s, B_d and B_f vanish, and energy is conserved.

3 SPATIAL DISCRETIZATION

Our aim is to discretize eqs (1) and (3) in a provably stable and accurate way, with a semi-discrete estimate that mimics (7). To do this we will use finite difference approximations satisfying an SBP rule with boundary and interface conditions enforced weakly using the SAT method (Kreiss & Scherer 1974, 1977; Strand 1994; Mattsson & Nordström 2004). We begin by describing the 1-D SBP operators and then move on to the Kronecker product construction of the 2-D operators. We conclude the section by describing the full 2-D discretization of eqs (1) and (3) and then stating a previously derived stability result.

Consider the 1-D domain Ω_y = [0, L_y] discretized into an evenly spaced grid of N_y + 1 points. We let the grid points be y_j = j h_y for j = 0, 1, …, N_y with spacing h_y = L_y/N_y. We let grid function |$\mathbf {p}^{T}=[p_0,\dots ,p_{N_y}]$| be the interpolation of function p with p_j = p(y_j). The first derivative of p is approximated as

$$ \begin{eqnarray*} \displaystyle\frac{\partial p}{ \partial y} &\approx \mathbf {D}_{1}\mathbf {p} = \mathbf {H}^{-1} \mathbf {Q} \mathbf {p}. \end{eqnarray*}$$

(8)

Here the finite difference matrix |$\mathbf {D}_{1}$| is of size (N_y + 1) × (N_y + 1). The matrix |$\mathbf {H}$| is diagonal and positive and can be thought of as a numerical quadrature rule (Hicken & Zingg 2013), namely,

$$ \begin{eqnarray*} \int _{0}^{L_{y}} p q\,\,\mathrm{ d}y \approx \mathbf {p}^{T} \mathbf {H} \mathbf {q}. \end{eqnarray*}$$

(9)

The matrix |$\mathbf {Q}$| is an almost skew symmetric matrix with the property that

$$ \begin{eqnarray*} \mathbf {Q} + \mathbf {Q}^{T} = \mathbf {B} = \mbox{diag}[-1,\,\,0,\,\,0\,\,\dots ,\,\,0,\,\,1]. \end{eqnarray*}$$

(10)

Operators with this structure are called SBP due to the identity

$$ \begin{eqnarray*} \mathbf {q}^{T} \mathbf {H} \mathbf {D}_{1} \mathbf {p} &=& \mathbf {q}^{T} \mathbf {Q} \mathbf {p} = \mathbf {q}^{T} \left(\mathbf {B} - \mathbf {Q}^{T}\right) \mathbf {p}\nonumber \\ & =& q_{N_{y}}p_{N_{y}} - q_{0}p_{0} - \mathbf {q}^{T} \mathbf {D}_{1}^{T} \mathbf {H} \mathbf {p}, \end{eqnarray*}$$

(11)

which discretely mimics integration by parts,

$$ \begin{eqnarray*} \int _{0}^{L_{y}} q \displaystyle\frac{\partial p}{ \partial y}\,\,dy &=& q\left(L_{y}\right)p\left(L_{y}\right) - q\left(0\right)p\left(0\right) \nonumber \\ && - \int _{0}^{L_{y}} \displaystyle\frac{\partial q}{ \partial y}p\,\,\mathrm{ d}y. \end{eqnarray*}$$

(12)

One approach to defining a second derivative operator is to apply |$\mathbf {D}_{1}$| twice:

$$ \begin{eqnarray*} \displaystyle\frac{\partial ^{2} p}{ \partial y^{2}} &\approx \mathbf {D}_{1}\mathbf {D}_{1}\mathbf {p}. \end{eqnarray*}$$

(13)

One downside of this is that it increases the bandwidth of the operator. Thus we instead prefer to use the compact SBP second derivative operators of Mattsson & Nordström (2004):

$$ \begin{eqnarray*} \displaystyle\frac{\partial ^{2} p}{ \partial y^{2}} &\approx \mathbf {D}_{2}\mathbf {p} = \mathbf {H}^{-1} \left(-\mathbf {M}+\mathbf {B} \mathbf {S}\right) \mathbf {p}. \end{eqnarray*}$$

(14)

The matrix |$\mathbf {M}$| is symmetric positive definite and can be thought of as approximating the inner product of derivatives:

$$ \begin{eqnarray*} \int _{0}^{L_{y}} \displaystyle\frac{\partial p}{\partial y} \displaystyle\frac{\partial q}{\partial y}\,\,\mathrm{ d}y \approx \mathbf {p}^{T} \mathbf {M} \mathbf {q}. \end{eqnarray*}$$

(15)

Matrix |$\mathbf {B}$| is as defined above and |$\mathbf {S}$| is an approximation of the first derivative; note that in general |$\mathbf {S} \ne \mathbf {D}_{1}$|⁠. We assume that |$\mathbf {H}$| in eqs (8) and (14) are the same, namely, the operators are compatible. Operator |$\mathbf {D}_{2}$| is called SBP since

$$ \begin{eqnarray*} \mathbf {p}^{T}\mathbf {H} \mathbf {D}_{2} \mathbf {q} = \mathbf {p}_{N_{y}}\left(\mathbf {S}\right)_{N_{y}} - \mathbf {p}_{0}\left(\mathbf {S}\right)_{0} -\mathbf {p}^{T}\mathbf {M} \mathbf {q} \end{eqnarray*}$$

(16)

discretely mimics the continuous identity

$$ \begin{eqnarray*} \int _{0}^{L_{y}} p \displaystyle\frac{\partial ^{2} q}{\partial y^{2}}\,\,\mathrm{ d}y &=& p\left(L_{y}\right) {\left.\displaystyle\frac{\partial q}{\partial y}\right|}_{y = L_{y}} - p\left(0\right) {\left.\displaystyle\frac{\partial q}{\partial y}\right|}_{y = 0} \nonumber \\ && - \int _{0}^{L_{y}} \displaystyle\frac{\partial p}{\partial y} \displaystyle\frac{\partial q}{\partial y}\,\,\mathrm{ d}y. \end{eqnarray*}$$

(17)

In this work we will exclusively consider the second-order accurate SBP operators, which are central difference operators in the interior and one-sided at the boundary. Note that the operators are second-order accurate in the interior but only first-order accurate at the boundary; however the global accuracy of these operator is 2 (Gustafsson 1975; Strand 1994; Mattsson & Nordström 2004). The operators we use are

$$ \begin{eqnarray*} \mathbf{D}_{1} &=& \displaystyle\frac{1}{h_{y}} \left[\begin{array}{cccccc} -1 & 1 & & & & \\ -\displaystyle\frac{1}{2} & 0 & \displaystyle\frac{1}{2} & & & \\ & -\displaystyle\frac{1}{2} & 0 & \displaystyle\frac{1}{2} & & \\ & & \ddots & \ddots & \ddots &\\ & & & -\displaystyle\frac{1}{2} & 0 & \displaystyle\frac{1}{2}\\ & & & & -1 & 1 \\ \end{array}\right], \nonumber \\ \mathbf {D}_{2} &=& \displaystyle\frac{1}{h_{y}^{2}} \left[\begin{array}{cccccc} 1 & -2 & 1 & & &\\ 1 & -2 & 1 & & &\\ & 1 & -2 & 1 & &\\ & & \ddots & \ddots & \ddots &\\ & & & 1 & -2 & 1\\ & & & 1 & -2 & 1 \end{array}\right], \end{eqnarray*}$$

(18)

where the SBP factors of the operators are

$$ \begin{eqnarray*} \mathbf {H} = \mbox{diag}\left[\displaystyle\frac{1}{2},\,\,1,\,\,1\,\,\dots ,\,\,1,\,\,\displaystyle\frac{1}{2}\right], \mathbf {Q} = {\begin{bmatrix}-\displaystyle\frac{1}{2} & \displaystyle\frac{1}{2}\\ -\displaystyle\frac{1}{2} & 0 & \displaystyle\frac{1}{2}\\ & \ddots & \ddots & \ddots \\ & & -\displaystyle\frac{1}{2} & 0 & \displaystyle\frac{1}{2} \\ & & & -\displaystyle\frac{1}{2} & \displaystyle\frac{1}{2} \end{bmatrix}},\nonumber \\ \end{eqnarray*}$$

(19)

$$ \begin{eqnarray*} \mathbf {S} = { \begin{bmatrix}-\displaystyle\frac{3}{2} & 2 & -\displaystyle\frac{1}{2} & \\ & 1 & \\ & & \ddots & \\ & & &1 \\ & & & \displaystyle\frac{1}{2} -2 & \displaystyle\frac{3}{2} \end{bmatrix}}, \mathbf {M} = {\begin{bmatrix}1 & -1 \\ -1 & 2 & -1\\ & \ddots & \ddots & \ddots \\ & & -1 & 2 & -1\\ & & & -1 & 1 \end{bmatrix}}. \end{eqnarray*}$$

(20)

The 1-D operators can be extended to multiple dimensions via the Kronecker product. The Kronecker product of matrices |$\mathbf {A}$| and |$\mathbf {C}$| is defined as

$$ \begin{equation*} \mathbf {A} \otimes \mathbf {C} = {\begin{bmatrix}a_{0,0}\mathbf {C} & \dots & a_{0,m}\mathbf {C}\\ \vdots & \ddots & \vdots \\ a_{m,0} \mathbf {C} & \dots & a_{m,n}\mathbf {C} \end{bmatrix}} \end{equation*}$$

where |$\mathbf {A}$| is of size m × n, |$\mathbf {C}$| is of size l × k and |$\mathbf {A} \otimes \mathbf {C}$| is of size ml × nk. We define the grid function of p as

$$ \begin{eqnarray*} {\bf \bar{p}} = [{\bf p}_0^T, {\bf p}_1^T, \, ...\, {\bf p}_{N_y}^T]^T \end{eqnarray*}$$

(21)

with |${\bf p}_i = [p_{i,0}, p_{i,1},\, ... \, p_{i,N_z}]^T$| for i = 0, ...N_y and p_{i, j} ≈ p(y_i, z_j). The derivative approximations are then:

$$ \begin{eqnarray*} \displaystyle\frac{\partial p}{\partial y} &\approx & \left({\bf D}_{1}^{(y)} \otimes {\bf I}^{(z)}\right){\bf \bar{p}} = \mathbf {\bar{D}}_{1}^{(y)} {\bf \bar{p}}, \nonumber \\ \displaystyle\frac{\partial p}{\partial z} &\approx & \left({\bf I}^{(y)} \otimes {\bf D}_{1}^{(z)}\right){\bf \bar{p}} = \mathbf {\bar{D}}_{1}^{(z)} {\bf \bar{p}}, \end{eqnarray*}$$

(22a)

$$ \begin{eqnarray*} \displaystyle\frac{\partial ^{2} p}{\partial y^{2}} &\approx & \left({\bf D}^{(y)}_{2} \otimes {\bf I}^{(z)}\right){\bf \bar{p}} = \mathbf {\bar{D}}_{2}^{(y)} {\bf \bar{p}}, \nonumber \\ \displaystyle\frac{\partial ^{2} p}{\partial z^{2}} &\approx & \left({\bf I}^{(y)} \otimes {\bf D}_{2}^{(z)}\right){\bf \bar{p}} = \mathbf {\bar{D}}_{2}^{(z)} {\bf \bar{p}}. \end{eqnarray*}$$

(22b)

Here |$\mathbf {I}^{(y)}$| and |$\mathbf {I}^{(z)}$| are identity matrices of size (N_y + 1) × (N_y + 1) and (N_z + 1) × (N_z + 1); the superscripts (y) and (z) in the derivative matrix indicate whether the operator is for the y or z dimensions, respectively.

With the above notation in place, we can now define the semi-discrete version of eqs (1) and (3) as (Virta & Mattsson 2014)

$$ \begin{eqnarray*} \rho \displaystyle\frac{\mathrm{ d}^2 \mathbf {\bar{u}}}{\mathrm{ d}t^2} &=& \mu _{1} \mathbf {\bar{D}}_{2}^{(y)} \mathbf {\bar{u}} + \mu _{2} \mathbf {\bar{D}}_{1}^{(y)} \mathbf {\bar{D}}_{1}^{(z)} \mathbf {\bar{u}} + \mu _{2} \mathbf {\bar{D}}_{1}^{(z)} \mathbf {\bar{D}}_{1}^{(y)} \mathbf {\bar{u}} + \mu _{3} \mathbf {\bar{D}}_{2}^{(z)}\mathbf {\bar{u}} + \mathbf {\bar{p}}_L \nonumber \\ && + \mathbf {\bar{p}}_R + \mathbf {\bar{p}}_B + \mathbf {\bar{p}}_T. \end{eqnarray*}$$

(23)

Here, the vectors |$\mathbf {\bar{p}}_L$|⁠, |$\mathbf {\bar{p}}_R$|⁠, |$\mathbf {\bar{p}}_B$|⁠, and |$\mathbf {\bar{p}}_T$| are penalty vectors that enforce the boundary and interface conditions. These vectors are defined as

$$ \begin{eqnarray*} \left(\mathbf {H}^{(y)} \otimes \mathbf {H}^{(z)}\right) \mathbf {\bar{p}}_{L} &=& \Big( \alpha \left(\mathbf {I}^{(y)} \otimes \mathbf {H}^{(z)}\right) - \left(\mu _{1} \mathbf {\bar{S}}^{(y)} + \mu _{2} \mathbf {\bar{D}}^{(z)}\right)^{T} \nonumber \\ && \times {\left(\mathbf {I}^{(y)} \otimes \mathbf {H}^{(z)}\right)} \Big) \mathbf {\bar{E}}_{L} \left(\mathbf {\bar{u}} - \mathbf {\bar{g}}_{L}\right), \end{eqnarray*}$$

(24a)

$$ \begin{eqnarray*} \left(\mathbf {H}^{(y)} \otimes \mathbf {H}^{(z)}\right) \mathbf {\bar{p}}_{R} &=& \Big( \alpha \left(\mathbf {I}^{(y)} \otimes \mathbf {H}^{(z)}\right) + \left(\mu _{1} \mathbf {\bar{S}}^{(y)} + \mu _{2} \mathbf {\bar{D}}^{(z)}\right)^{T}\nonumber \\ && \times {\left(\mathbf {I}^{(y)} \otimes \mathbf {H}^{(z)}\right)} \Big) \mathbf {\bar{E}}_{R} \left(\mathbf {\bar{u}} - \mathbf {\bar{g}}_{R}\right), \end{eqnarray*}$$

(24b)

$$ \begin{eqnarray*} \left(\mathbf {H}^{(y)} \otimes \mathbf {H}^{(z)}\right) \mathbf {\bar{p}}_{B} &= {\left(\mathbf {H}^{(y)} \otimes \mathbf {I}^{(z)}\right)} \mathbf {\bar{E}}_{B} \left(\mu _{2} \mathbf {\bar{D}}^{(y)} + \mu _{3} \mathbf {\bar{S}}^{(z)}\right) \mathbf {\bar{u}}, \end{eqnarray*}$$

(24c)

$$ \begin{eqnarray*} \left(\mathbf {H}^{(y)} \otimes \mathbf {H}^{(z)}\right) \mathbf {\bar{p}}_{T} &= - {\left(\mathbf {H}^{(y)} \otimes \mathbf {I}^{(z)}\right)}\mathbf {\bar{E}}_{T} \left(\mu _{2} \mathbf {\bar{D}}^{(y)} + \mu _{3} \mathbf {\bar{S}}^{(z)}\right) \mathbf {\bar{u}}. \end{eqnarray*}$$

(24d)

Here the vectors |$\mathbf {\bar{g}}_{L}$| and |$\mathbf {\bar{g}}_{R}$| are the grid functions which are zero everywhere except for along the left and right boundaries where they take the values of g_L and g_R, respectively (see eq. 3). The matrices |$\mathbf {\bar{E}}_{L}$|⁠, |$\mathbf {\bar{E}}_{R}$|⁠, |$\mathbf {\bar{E}}_{B}$| and |$\mathbf {\bar{E}}_{T}$| zero out all values in a vector except those along the left, right, bottom and top boundaries, respectively, and are defined as

$$ \begin{eqnarray*} \mathbf {\bar{E}}_{L} &=& \mbox{diag}(1,~0,~\dots ,~0) \otimes \mathbf {I}^{(z)}, \mathbf {\bar{E}}_{R} \nonumber \\ & = & \mbox{diag}(0,~\dots ,~0,~1) \otimes \mathbf {I}^{(z)}, \end{eqnarray*}$$

(25a)

$$ \begin{eqnarray*} \mathbf {\bar{E}}_{B} &=& \mathbf {I}^{(y)} \otimes \mbox{diag}(0,~\dots ,~0,~1), \mathbf {\bar{E}}_{T} \nonumber \\ &=& \mathbf {I}^{(y)} \otimes \mbox{diag}(1,~0,~\dots ,~0). \end{eqnarray*}$$

(25b)

The 2-D boundary derivative matrices are

$$ \begin{eqnarray*} \mathbf {\bar{S}}^{(y)} = \mathbf {S}^{(y)} \otimes \mathbf {I}^{(z)}, \mathbf {\bar{S}}^{(z)} = \mathbf {S}^{(z)} \otimes \mathbf {I}^{(z)}. \end{eqnarray*}$$

(26a)

Penalty terms |$\mathbf {\bar{p}}_{B}$| and |$\mathbf {\bar{p}}_{T}$| enforce the free surface boundary conditions whereas |$\mathbf {\bar{p}}_{L}$| and |$\mathbf {\bar{p}}_{R}$| enforce the Dirchlet boundary and fault interface conditions (hence the need to subtract off data from the solution vector). The scalar parameter α in |$\mathbf {\bar{p}}_{L}$| and |$\mathbf {\bar{p}}_{R}$| needs to be sufficiently large (in magnitude) so that the discretization is energy stable. For the second-order accurate operators used here, the results of Virta & Mattsson (2014) (reduced to the case of constant coefficients) show that α must satisfy

$$ \begin{eqnarray*} &&\alpha \le -\displaystyle\frac{99}{36}\displaystyle\frac{\mu _{1}^{2}}{ \lambda h_y} -\displaystyle\frac{2 \mu _{3}^{2}}{\lambda h_y},\nonumber \\ &&{\lambda} = \tfrac{1}{2}(\mu _1+\mu _3) -\sqrt{(\mu _1-\mu _3)^2 + 4\mu _{2}^{2}}. \end{eqnarray*}$$

(27)

With zero boundary data, g_L = g_R = 0, Virta & Mattsson (2014) derive an energy estimate for the numerical solution to the semi-discrete equations, showing the scheme is energy stable.

4 FRICTIONAL FRAMEWORK

The displacements and tractions on the two sides of a fault interface, located at y = 0 in our model, are related to one another via a nonlinear friction law that enforces continuity of traction while allowing for a jump in displacement. We define the slip velocity or the time derivative of the jump in displacement across the fault by

$$ \begin{eqnarray*} V(z, t) &= \displaystyle\frac{\partial \Delta u(z,t)}{\partial t}, \end{eqnarray*}$$

(28a)

$$ \begin{eqnarray*} \Delta u(z,t) &= \lim _{\epsilon \rightarrow 0^{+}} \left(u(\epsilon , z, t) - u(-\epsilon , z, t)\right), \end{eqnarray*}$$

(28b)

and the shear stress on the fault by

$$ \begin{eqnarray*} \tau = \sigma _{xy}(0,z,t) = \left(\mu _{1} \displaystyle\frac{\partial u}{\partial y} + \mu _{2} \displaystyle\frac{\partial u}{\partial z}\right)\Big \vert _{y = 0}, \end{eqnarray*}$$

(29)

that is, the component of traction in the x-direction, on the y ≤ 0 side of the interface, that comes from the y ≥ 0 side. Rate-and-state friction relates the shear stress τ on the fault to a nonlinear function of the slip velocity V and a state variable Ψ which obeys a local ordinary differential equation that tracks the history of sliding (Dieterich 1979; Marone 1998):

$$ \begin{eqnarray*} \tau &= F(V, \Psi ), \end{eqnarray*}$$

(30a)

$$ \begin{eqnarray*} \displaystyle\frac{\mathrm{ d} \Psi }{\mathrm{ d} t} &= G(V, \Psi ). \end{eqnarray*}$$

(30b)

These relationships along with continuity of traction, that is, Δσ_xy = 0 across the fault, fully specify the problem. The specific forms of F and G we use are as follows:

$$ \begin{eqnarray*} F(V, \Psi ) &= a \sigma _{n} \sinh ^{-1}\left(\displaystyle\frac{V}{2V_{0}} e^{\displaystyle\frac{\Psi }{a}}\right), \end{eqnarray*}$$

(31a)

$$ \begin{eqnarray*} G(V,\Psi ) &= \displaystyle\frac{b V_{0}}{D_{c}} \left(e^{\displaystyle\frac{f_{0} - \Psi }{b}} - \displaystyle\frac{V}{V_{0}}\right), \end{eqnarray*}$$

(31b)

where f₀ is a reference friction coefficient for steady sliding at slip velocity V₀, a and b are dimensionless parameters characterizing the direct and state evolution effects, respectively, σ_n is the effective normal stress on the fault and D_c is the state evolution distance.

An important feature of the friction law is that even though the governing equations are linear in the volume, friction law (31a) is nonlinear. This nonlinearity poses no computational challenge if explicit time integration is used for semi-discretization (23), as the friction law only enters on the right-hand side of the equation; see for instance Kozdon et al. (2012). The problem with using explicit time integration for earthquake cycle simulations is that the CFL restriction will lead to a very small time step (on the order of milliseconds with realistic material parameters) which would make long time simulations (hundreds of years) impractical. One approach would be to use implicit time stepping when the slip velocity V along the whole fault is low to thus “step over” the extremely low frequency waves. The problem with this approach is that the friction law (31a) then leads to a large nonlinear system of equations that must be solved. Thus here, following Erickson & Dunham (2014), we set the inertial term |$\mathrm{ d}^2 \mathbf {\bar{u}}/\mathrm{ d}t^2$| in the semi-discretization (23) to zero. With this semi-discretization (23) then becomes a linear system of the form

$$ \begin{eqnarray*} \mathbf {\bar{A}} \mathbf {\bar{u}} = \mathbf {\bar{b}}(\mathbf {\Delta u}, t). \end{eqnarray*}$$

(32)

Here |$\mathbf {\bar{A}}$| is a matrix of size N_p × N_p and |$\mathbf {\bar{b}}(\mathbf {\Delta u}, t)$| is a vector of size N_p where in both cases N_p = (N_y + 1)(N_z + 1) is the total number of grid points. The vector |$\mathbf {\bar{b}}(\mathbf {\Delta u}, t)$| incorporates the boundary conditions, which due to the friction law and outer boundary depend on both t and Δu. Note that in semi-discretization (23), symmetry implies that g_L = Δu/2, namely, the jump in displacement on the fault is accommodated equally on both sides.

By zeroing out the inertial terms we are then saying that changes in displacement on the fault (and outer boundary) instantaneously modify displacements in the interior. This assumption is valid when the magnitude of the slip velocity is low, |V| ≪ 1, but for higher sliding velocities waves must be approximated in some way for the problem to remain relevant and well-posed. Here we use the radiation damping approximation (Rice 1993). In this approach waves that result from slip on the fault are assumed to emit shear waves that propagate normal to the fault. The effect of this is that shear stress on the fault is decreased by a factor of ηV where |$\eta = \sqrt{\mu _{1} \rho } / 2$| is half the shear-wave impedance. With this, the friction law is modified to

$$ \begin{eqnarray*} \tau _{\mathrm{ qs}} - \eta V = F(V, \Psi ), \end{eqnarray*}$$

(33)

where τ_qs is the “quasi-static” shear stress (computed via eq. 29), based on solving (23) without inertial terms.

In this formulation, time enters the equation through the state evolution eq. (30b) and when Δu is updated using (28a). Given a value of Δu and Ψ all that remains to be determined is V (since G can be evaluated once V and Ψ are known). To determine |$\mathbf {V}$| the following approach is used at a time t given |$\mathbf {\Delta u}$| and |$\boldsymbol{\Psi }$| (here we use vector notation to denote that these quantities are grid function along the fault).

1. The linear system (32) is solved for |$\mathbf {\bar{u}}$|

2. The displacement vector |$\mathbf {\bar{u}}$| is then used to compute |$\mathbf {\tau }$| as

   $$ \begin{eqnarray*} \boldsymbol{\tau }_{\mathrm{ qs}} = \left([1\,\,0\,\,\dots \,\,0] \otimes \mathbf {I}^{(z)}\right) \left(\mu _{1} \mathbf {\bar{S}}^{(y)} \mathbf {\bar{u}} + \mu _{2} \mathbf {\bar{D}}^{(z)}\right) \end{eqnarray*}$$

   (34)
3. At each grid point along the fault the nonlinear system

   $$ \begin{eqnarray*} {\left[\tau _{\mathrm{ qs}}\right]}_{i} - \eta V_{i} = F\left(V_{i}, \Psi _{i}\right), \end{eqnarray*}$$

   (35)

   is solved for V_i. Here V_i, [τ_qs]_i, and Ψ_i are the values of these variables at each of the grids points with i = 0, 1, …, N_z.

The ODEs are then integrated in time using an adaptive time step Runge–Kutta method.

4.1 Adaptive time stepping

As mentioned above, explicit time integration requires time steps on the order of milliseconds due to restrictions of the CFL condition. To illustrate this point consider the simple example of a 1-D explicit time-marching discretization of the transport equation, in this setting the condition is

$$ \begin{align*} \left|\dfrac{\alpha \Delta t}{\Delta x} \right|\le 1, \end{align*}$$

Where α is the propagation velocity of the solution wave, Δt is the time step size and Δx is the spatial resolution. For seismic waves, propagation velocities tend to range between 2 and 8 km s⁻¹ (Shearer 2009). Take α = 8 km s⁻¹ and grid spacing Δx = 0.1 km, for example,

$$ \begin{align*} && \text{ The CFL condition requires } \quad \left|\dfrac{8 \Delta t}{0.1} \right|\le 1,\nonumber \\ && \text{ leading to a time step }\Delta t \le 0.0125 \text{ s.} \end{align*}$$

Simulating thousands of years at a step size on the order of milliseconds is computationally unfeasible. In Virta & Mattsson (2014) the anisotropic acoustic wave equation is modelled within heterogeneous media with an explicit time-integration scheme. The time step size is taken to be |$\dfrac{1}{2}\dfrac{\min (\Delta x,\Delta y, \Delta z)}{1700}$|⁠, where Δx, Δy and Δz are the grid spacings in the x, y and z directions. In the context of the grid spacings used in this paper this would correspond to a time step on the order of 10⁻⁵. In Komatitsch & Vilotte (1998) where elastic-wave propagation is modelled with the spectral element method, the time steps range from 0.25 to 0.5 ms.

In order to quantify the timescales that the adaptive time stepping is able to capture, we calculated the time step sizes for all of the simulations presented in Section 7, after these were out of the spin-up period. The smallest step size amidst all of the simulations was approximately 0.275 milliseconds, while the largest step size was approximately 6 yr. The difference between these is over 10 orders of magnitude. In general, smallest step sizes ranged between approximately 0.275 milliseconds and 0.68 milliseconds and the largest ranged between approximately 2 and 6 years. Comparing the step sizes employed in (Komatitsch & Vilotte 1998) and (Virta & Mattsson 2014) to those in this paper, we see that while the smallest step sizes that occur are in line with those used in explicit time-integration schemes, the largest are over 10 orders of magnitude larger.

5 DETERMINATION OF MAXIMUM STABLE GRID SPACING: LINEAR STABILITY ANALYSIS OF FRICTIONAL SLIDING

One of the important considerations in fracture modelling, such as earthquake cycle simulations, is the determination of the maximum stable grid spacing required (along the fault) so that numerical errors do not trigger unstable slip. That is, a well-posed fracture mechanics problem will have some critical wavelength h* such that perturbations which have wavelengths h* will decay in time whereas perturbations with wavelength greater than h* will grow (leading to unstable sliding). The implication for grid spacing in a finite difference method is that h* must be resolved with at least a few grid points so that the smallest wavelength discrete solutions decay; there is an additional length scale, known as cohesive, or process zone size, that must also be resolved and this is discussed at the end of this section.

In order to determine h* we extend the linear stability analysis of Ranjith & Gao (2007) to the anisotropic case. We considering antiplane sliding of two identical anisotropic elastic half-spaces separated by a frictional fault at y = 0.

We Laplace transform the equilibrium version of eq. (1) in time to obtain

$$ \begin{eqnarray*} 0=\mu _1 \dfrac{\partial ^2 \hat{u}}{\partial y^2} + 2\mu _2 \dfrac{\partial ^2 \hat{u}}{\partial y \partial z} +\mu _3 \dfrac{\partial ^2 \hat{u}}{\partial z^2}. \end{eqnarray*}$$

(36)

Letting the solution to eq. (36) be of the form

$$ \begin{eqnarray*} \hat{u}(y,z,p) = \hat{\mathcal {U}}(y,k,p)\mathrm{ e}^{ikz}, \end{eqnarray*}$$

(37)

where p is the Laplace transformed variable, we then get

$$ \begin{eqnarray*} 0 = \mu _1 \dfrac{\partial ^2 \hat{\mathcal {U}}}{\partial y^2} + 2 ik \mu _2 \dfrac{\partial \hat{\mathcal {U}}}{\partial y} +(ik)^2 \mu _3 \hat{\mathcal {U}}. \end{eqnarray*}$$

(38)

Solutions of the ordinary differential equation (38) are of the form

$$ \begin{eqnarray*} \hat{\mathcal {U}}(y,k,p) = \left\lbrace \begin{array}{@{}l@{\quad }l@{}}\hat{\mathcal {U}}^{(+)}(k,p)\mathrm{ e}^{\alpha ^{(+)}y}, \quad y\gt 0,\\ \hat{\mathcal {U}}^{(-)}(k,p)\mathrm{ e}^{\alpha ^{(-)}y}, \quad y\lt 0. \end{array}\right. \end{eqnarray*}$$

(39)

Here we have used the superscript (+) to denote the positive side of the fault (y > 0), and (−) to denote the side y < 0. The characteristic roots α^(±) are found by substituting in solution (39) into ordinary differential equation (38) and solving the resulting quadratic equation:

$$ \begin{eqnarray*} \alpha ^{(\pm )}=\displaystyle\frac{-ik\mu _2\mp \sqrt{k^2(\mu ^*)^2}}{\mu _1}, \quad \mu ^*=\sqrt{\det {{\begin{bmatrix}\mu _1 &\mu _2\\ \mu _2 & \mu _3 \end{bmatrix}}}}; \end{eqnarray*}$$

(40)

here we have chosen the root on each side of the fault so that α^(±)y has a negative real part and the thus the solution decays as |y| → ∞. Putting this all together then yields

$$ \begin{eqnarray*} \hat{u}(y,z,p)= \left\lbrace \begin{array}{@{}l@{\quad }l@{}}\hat{\mathcal {U}}^{(+)}(k,p)\mathrm{ e}^{\alpha ^{(+)}y+ikz}, \quad y\gt 0,\\ \hat{\mathcal {U}}^{(+)}(k,p)\mathrm{ e}^{\alpha ^{(-)}y+ikz}, \quad y\lt 0, \end{array}\right. \end{eqnarray*}$$

(41)

or upon transforming back from Laplace space

$$ \begin{eqnarray*} u(y,z,t)= \left\lbrace \begin{array}{@{}l@{\quad }l@{}}\mathcal {U}^{(+)}(k,t)e^{\alpha ^{(+)}y+ikz}, \quad y\gt 0,\\ \mathcal {U}^{(+)}(k,t)e^{\alpha ^{(-)}y+ikz}, \quad y\lt 0. \end{array}\right. \end{eqnarray*}$$

(42)

The Laplace transform of traction on the two sides of the fault at y = 0 is

$$ \begin{eqnarray*} \hat{\sigma }_{xz}^{(\pm )} = \hat{T}^{(\pm )} \mathrm{ e}^{ikz}, \quad \hat{T}^{(\pm )} = \left(\mu _1 \alpha ^{(\pm )} + ik\mu _2\right)\hat{\mathcal {U}^{(\pm )}}(k,t). \end{eqnarray*}$$

(43)

Continuity of traction implies that |$\hat{T}^{(+)} = \hat{T}^{(-)}$|⁠, which after some simplification gives

$$ \begin{eqnarray*} -|k| \mu ^{*} \hat{\mathcal {U}}^{(+)}(k, p) = |k| \mu ^{*} \hat{\mathcal {U}}^{(-)}(k, p). \end{eqnarray*}$$

(44)

The jump in displacement across the fault is

$$ \begin{eqnarray*} && \hat{u}(0^{+},z,p) - \hat{u}(0^{-},z,p) = \hat{D}(k, p) \mathrm{ e}^{ikz}, \nonumber \\ && \hat{D}(k, p) = \hat{\mathcal {U}}^{(+)}(k,p) - \hat{\mathcal {U}}^{(-)}(k,p). \end{eqnarray*}$$

(45)

Using |$\hat{D}$| along with continuity of traction allows us then to relate slip to traction:

$$ \begin{eqnarray*} \hat{T}(k,p) = \displaystyle\frac{\hat{T}^{(+)}(k,p) + \hat{T}^{(-)}(k,p)}{2} = -\displaystyle\frac{|k|}{2}\mu ^{*} \hat{D}(k,p). \end{eqnarray*}$$

(46)

Laplace transforming of the time derivative of the linearized rate-and-state friction law is (Ranjith & Gao 2007)

$$ \begin{eqnarray*} \left(p+\displaystyle\frac{V_0}{D_c}\right)\hat{T}=\displaystyle\frac{\sigma _n}{V_0}\left(ap-(b-a)\displaystyle\frac{V_0}{D_{c}}\right)p\hat{D}. \end{eqnarray*}$$

(47)

Substitution of eq. (46) into eq. (47) yields the quadratic

$$ \begin{eqnarray*} \displaystyle\frac{\sigma _n}{V_0}ap^2 +p\left(\displaystyle\frac{|k|}{2}\mu ^* -\displaystyle\frac{(b-a)\sigma _n}{D_c}\right)+\displaystyle\frac{V_0}{D_c}\displaystyle\frac{|k|}{2}\mu ^* =0. \end{eqnarray*}$$

(48)

The system will undergo Hopf bifurcation when roots p cross the imaginary axis, which will occur when |k| is less than the critical wave number k_cr:

$$ \begin{eqnarray*} |k_{\mathrm{ cr}}|=\dfrac{2(b-a)\sigma _n}{D_c\mu ^*}. \end{eqnarray*}$$

(49)

In terms of wavelength, this corresponds to the critical wavelength

$$ \begin{eqnarray*} h^* &= \displaystyle\frac{2\pi }{|k_{\mathrm{ cr}}|} = \displaystyle\frac{\pi \mu ^*D_c}{(b-a)\sigma _n}. \end{eqnarray*}$$

(50)

This then implies that we want our grid spacing to be smaller than h* so that numerical noise does not trigger ruptures.

As noted above, Ampuero & Rubin (2008) derive an even smaller length scale called the cohesive zone, which, for the anisotropic problem, we interpret to be

$$ \begin{eqnarray*} L_b = \dfrac{\mu ^{*}D_c}{\sigma _n b}. \end{eqnarray*}$$

(51)

The cohesive zone length L_b corresponds to the spatial length scale over which the shear stress drops from its peak to residual values at the propagating rupture front; numerical studies of quasi-dynamic earthquake cycle simulations in isotropic materials observe that L_b must be resolved with at least one grid point (Ampuero & Rubin 2008). In our simulations we resolve L_b with at least 5 grid points. To ensure that this grid spacing is sufficient, we doubled the number of grid points so that L_b and h* were resolved with over 10 and 120 grid points respectively. Comparison to simulations with doubled resolution indicates that resolving L_b with over 5 grid points is adequate; see Appendix  A.

6 CONVERGENCE TESTS

We verify our numerical method via the method of manufactured solutions (Roache 1998). In this approach source terms are added to eqs (1) and (3) so that a known function can be used as an analytic solution. Namely, we let the exact displacement |$\hat{u}$| be given as:

$$ \begin{eqnarray*} \hat{u}(t,y,z)=\displaystyle\frac{\delta }{2}K(t)\phi (y,z)+\displaystyle\frac{V_{p}t}{2}\left(1-\phi (y,z)\right)+\displaystyle\frac{\tau ^{\infty }}{\mu _1}y, \end{eqnarray*}$$

(52)

where K(t) and ϕ(y, z) are functions which will determine the temporal and spatial dependence of the solution. These functions will be chosen so that the solution exhibits both an interseismic (slow) and a coseismic (fast) phase. Namely, there will be an initial interseismic phase, followed by a single coseismic phase, and another interseismic phase; this allows us to verify the ability of our time stepping method to integrate accurately through these different phases. Parameters V_p and τ^∞ are the plate rate and magnitude of remote stress and are taken to be constant; see Table 1. The parameter δ is the total slip during the coseismic phase, and we take it to be equal to δ = (V_p + V_min)t_e, where t_e is the time of the coseismic event and V_min is the minimum slip velocity; see Table 1.

The spatial dependency of the manufactured solution is given by

$$ \begin{eqnarray*} \phi (y,z) = \dfrac{H(H+y)}{(H+y)^2+z^2}, \end{eqnarray*}$$

(53)

where H is a locking depth given in Table 1. When evaluated along the fault (at y = 0) ϕ takes the form of a normalized Lorentzian distribution. The temporal dependency of the manufactured solution is

$$ \begin{eqnarray*} K(t)=\displaystyle\frac{1}{\pi }\left[\arctan \left(\displaystyle\frac{t-t_e}{t_w}\right)+\displaystyle\frac{\pi }{2}\right]+\displaystyle\frac{V_{\min }}{\delta }t, \end{eqnarray*}$$

(54)

and is designed to test the adaptive time-stepping of the numerical scheme. The system is first loaded at the plate rate of V_p, this corresponds to the interseismic period. When the system reaches event time t_e, slip velocity increases over many orders of magnitude, simulating a rupture event with duration t_w. Velocity then returns to its minimum rate, V_min, for the rest of the simulation.

Table 1.

Parameters used in manufactured solution convergence tests.

Parameter	Definition	Value
Ly	Fault domain length	72 km
Lz	Off-fault domain length	72 km
H	Locking depth	12 km
μ1	Material stiffness parameter	36 GPa
μ2	Material stiffness parameter	variable GPa
μ3	Material stiffness parameter	variable GPa
ρ	Density	2800 kg/m3
σn	Normal stress in fault	50 MPa
τ∞	Remote shear stress	31.73 MPa
tf	Final simulation time	70 yr
te	Event nucleation time	35 yr
tw	Timescale for event duration	10 s
a	Rate-and-state parameter	0.015
b	Rate-and-state parameter	0.02
Dc	critical slip distance	0.2 m
Vp	Plate rate	10−9 m s−1
Vmin	Minimum slip velocity	10−12 m s−1
V0	Reference velocity	10−6 m s−1
f0	Reference friction coefficient	0.6

The exact solution sets the initial data for differential eqs (28a) and (30b) and allows us to solve for the exact shear stress |$\hat{\tau }_{\mathrm{ qs}}$| via (29) and the exact slip velocity |$\hat{V}$| via (28). Plugging these into (33) allows us to solve for |$\hat{\psi }$|⁠, namely

$$ \begin{eqnarray*} \hat{\psi } = a\ln \left[ \displaystyle\frac{2V_0}{\hat{V}}\sinh \left(\dfrac{\hat{\tau }_{qs}-\eta \hat{V}}{\sigma _n a}\right)\right]. \end{eqnarray*}$$

(55)

The boundary data g_L(t, z) is obtained via integration of the ODEs, as detailed in Section 4. Note that as done in Erickson & Dunham (2014), we must add a source term to (30b), that is, to update state evolution we now numerically integrate

$$ \begin{eqnarray*} \displaystyle\frac{\partial \psi }{\partial t}= G(V,\psi )+ \displaystyle\frac{\partial \hat{\psi }}{\partial t}- G(\hat{V},\hat{\psi }). \end{eqnarray*}$$

(56)

The manufactured solution we have chosen does not satisfy the traction free boundary condition, and thus we instead enforce the top and bottom stress boundary conditions:

$$ \begin{eqnarray*} \sigma _{xz}(t,y,0) &=& \left[\mu _2 \displaystyle\frac{\partial \hat{u}}{\partial y} + \mu _3 \displaystyle\frac{\partial \hat{u}}{\partial z}\right] \bigg |_{z=0}, \nonumber \\ \sigma _{xz}(t,y,L_z) &=& \left[\mu _2 \displaystyle\frac{\partial \hat{u}}{\partial y} + \mu _3 \displaystyle\frac{\partial \hat{u}}{\partial z}\right]\bigg |_{z=L_z}. \end{eqnarray*}$$

(57a)

Similarly, the remote boundary data is defined by the manufactured solution,

$$ \begin{eqnarray*} g_R(t,z) = \hat{u}(t,L_y,z). \end{eqnarray*}$$

(58)

Because our main focus is to explore the effects of anisotropy within a homogeneous medium, we run convergence tests for both the orthotropic (μ₂ = 0) and fully anisotropic (μ₂ ≠ 0) cases with constant coefficients, and verify that the numerical solution converges to the exact solution at the expected rate for a second-order accurate method. At the end of each simulation, we compute the relative error in the discrete H-norm, given by

$$ \begin{eqnarray*} \text{Error}_H(h)=\dfrac{\left\Vert \mathbf {\hat{u}}-\mathbf {u}\right\Vert _H}{\left\Vert \mathbf {\hat{u}}\right\Vert _H}. \end{eqnarray*}$$

(59)

All the parameter values used in the convergence test simulations are located in Table 1. Tables 2 and 3 show the successive relative errors and convergence rates under mesh refinement for the orthotropic and fully anisotropic cases respectively.

Table 2.

Relative error for the orthotropic case (μ₃ = 24 GPa, μ₂ = 0 GPa), computed in the discrete |$\mathbf {H}$| norm with N = N_y = N_z. The rate of convergence approaches 2 under mesh refinement.

N	Error (h)	Rate
24 + 1	2.2541 × 10−2	–
25 + 1	6.0595 × 10−3	1.89527880
26 + 1	1.5770 × 10−3	1.94205097
27 + 1	4.0235 × 10−4	1.97063645
28 + 1	1.0072 × 10−4	1.99809893

Table 3.

Relative error for the fully anisotropic case (μ₂ = 18 GPa, μ₃ = 36 GPa), computed in the discrete |$\mathbf {H}$| norm with N = N_y = N_z. The rate of convergence approaches 2 under mesh refinement.

N	Error (h)	Rate
24 + 1	3.3244 × 10−2	–
25 + 1	8.9966 × 10−3	1.88564247
26 + 1	2.3099 × 10−3	1.96151856
27 + 1	5.82 × 10−4	1.98782444

Table 4.

Parameters used in the parameter-varying study.

Parameter	Definition	Value
Ly	Fault domain length	72 km
Lz	Off-fault domain length	72 or 120 km (see Appendix  B)
μ1	Material stiffness parameter	variable GPa
μ2	Material stiffness parameter	variable GPa
μ3	Material stiffness parameter	variable GPa
ρ	Density	2800 kg m−3
σn	Normal stress in fault	50 MPa
τ∞	Remote shear stress	31.73 MPa
a	Rate-and-state parameter	0.015
b	Rate-and-state parameter	depth variable
Dc	Critical slip distance	8 mm
Vp	Plate rate	10−9 m s−1
V0	Reference velocity	10−6 m s−1
f0	Reference friction coefficient	0.6

7 RESULTS OF PARAMETER-VARYING STUDY

Here we use the numerical scheme detailed in Sections 3 and 4 to study earthquake cycles in anisotropic media. We begin by considering an orthotropic medium (μ₂ = 0) and conduct a parameter study by varying μ₁ and μ₃ (holding all other parameters fixed). Rate-and-state friction parameters a and b vary with depth, as illustrated in Fig. 1. Negative values of a − b correspond to the velocity-weakening (seismogenic) zone where earthquakes occur. At depths for which a − b is positive, below ∼12 km depth, the fault undergoes steady-state sliding. We find that anisotropy influences the periodicity of a simulation, the length of the interseismic period, and that many of the simulations host aseismic transient events. In Section 7.3 we investigate the relationship between these transients, nucleation zone and interseismic creep. Due to the large scope of the parameter study, and the complexity of the results, we begin by studying the orthotropic problem and examine some illustrative examples. With the groundwork in place to better understand the entire set of results, we discuss all of our findings in Section 7.4. The entire parameter study is summarized in Table 5. In Section 7.5 we present preliminary results from the study of the fully anisotropic problem, and Section 7.6 illustrates the effects of anisotropy on surface velocity profiles.

Figure 1.

Depth variation of frictional parameters a and b shown down to a depth of 24 km. Below this depth the parameters remain constant, namely a − b = a = 0.015.

Table 5.

μ₁ values are located in the first column, while μ₃ values are in the first row. The interior cells contain the μ* value that corresponds to each μ₁, μ₃ pair. Colour indicates the period of a simulation, yellow is period one, red is period two and blue is period three. Bold font denotes a simulation that hosts aseismic transients. Circled cells are simulations for which a period change occurs for parameter-differences in simulations with the same μ* value.

Even though inertial effects are not considered in our simulations, it is useful to consider the wave speeds for an anisotropic medium in order to connect the model parameters with observations. Since the material stiffness matrix |$\mathbf {M}_{\mu }$| is symmetric, positive definite, it can be diagonalized as

$$ \begin{eqnarray*} \mathbf {M}_{\mu } &= {\begin{bmatrix}\mu _1 & \mu _2 \\ \mu _2 & \mu _3 \\ \end{bmatrix}} = \mathbf {V}{\begin{bmatrix}\lambda _1& 0 \\ 0& \lambda _2 \\ \end{bmatrix}}\mathbf {V}^T, \end{eqnarray*}$$

(60)

where |$\mathbf {V}$| has the orthogonal eigenvectors of |$\mathbf {M}_{\mu }$| as its columns and λ₁, λ₂ > 0. This means that the orthogonally split shear waves have fast and slow wave speeds

$$ \begin{eqnarray*} c_{\text{fast}}&=\sqrt{\displaystyle\frac{\max \lbrace \lambda _1,\lambda _{2}\rbrace }{\rho }}, \end{eqnarray*}$$

(61a)

$$ \begin{eqnarray*} c_{\text{slow}}&=\sqrt{\displaystyle\frac{\min \lbrace \lambda _1,\lambda _{2}\rbrace }{\rho }}. \end{eqnarray*}$$

(61b)

The fast wave travels in the direction of the eigenvector associated with the maximum eigenvalue and the slow wave travels in the direction of the eigenvector associated with the minimum eigenvalue.

Anisotropy measurements from shear wave splitting techniques used following the M7.1 Hector Mine earthquake (which occurred along a strike-slip fault) find anisotropy confined to the upper 2–3 km depth of the fault, with average fast directions oriented between fault parallel (within the horizontal plane) and parallel to the direction of regional maximum compressive stress (Cochran et al. 2003). Fast direction oriented with the fault-perpendicular direction, however, have also been observed (Stuart et al. 2002). Cochran et al. (2003) compute an apparent crack density ϵ = v_sδt/L, where v_s is the fast shear wave velocity, and δt/L is a path-normalized delay time between arrivals of the fast and slow shear waves, and report values of ϵ to be approximately 5 per cent, and generally less than 10 per cent, regardless of region. Thus the relationship between wave speeds is

$$ \begin{eqnarray*} c_{\mathrm{ slow}}(1 + \epsilon ) = c_{\mathrm{ fast}}, \end{eqnarray*}$$

(62)

which translates to a relationship between shear moduli given by

$$ \begin{eqnarray*} \max \lbrace \lambda _1,\lambda _{2}\rbrace / \min \lbrace \lambda _1,\lambda _{2}\rbrace = (1 + \epsilon )^2. \end{eqnarray*}$$

(63)

A maximum value of ϵ of 10 per cent corresponds to an approximate 20 per cent difference in shear moduli. We consider parameter values both within and outside this range in order to explore the full range of possible effects of anisotropy on the earthquake cycle.

7.1 Orthotropic Anisotropy

For orthotropic anisotropy (μ₂ = 0), two orthogonally split shear waves travel in either the fault-perpendicular (y-) or fault-parallel (z-) directions. If μ₁ > μ₃, the material is called horizontal transversely isotropic (HTI) and the fast wave travels in the fault perpendicular direction. When μ₁ < μ₃ the material is called vertical transversely isotropic (VTI) and the fast wave travels in the fault parallel direction.

For the orthotropic problem we first examined the effects of anisotropy by holding μ₁ constant and decreasing μ₃ (or vice versa). This corresponds to increasing the degree of HTI or VTI anisotropy. For a related isotropic problem, we consider several cases. One possibility is to choose the isotropic shear modulus, μ, to be either μ₁ or μ₃. An alternative isotropic reference case would be for a given to choose |$\mu = \mu ^{*} = \sqrt{\mu _1 \mu _3}$|⁠. In the text that follows, each of these choices will be used as relevant to the discussion.

7.2 Simulation Results

We first illustrate our findings in Fig. 2, where we show results for both an HTI and VTI simulation with μ* = 24 GPa, along with two isotropic reference cases. In these snapshots, slip is plotted over a sequence of earthquakes spanning about 1500 yr, where contours in blue represent the interseismic period, when the maximum slip rate (taken over the fault) is less than 1 mm s⁻¹. Slip is plotted in red contours every second during a quasi-dynamic event, when the maximum slip rate is greater than 1 mm s⁻¹. Fig. 2(a) is the HTI simulation, with μ₁ = 34 GPa and μ₃ = 16.95 GPa. Fig. 2(b) shows the VTI simulation, where μ₁ = 16.95 GPa and μ₃ = 34 GPa. Fig. 2(c) corresponds to an isotropic reference case, where μ₁ = μ₃ = 24 GPa (to ensure that μ* = 24 GPa as in each of the anisotropic simulations). Fig. 2(d) is the isotropic reference case for which μ₁ = μ₃ = 34 GPa. Comparing the snapshots of cumulative slip for the first three simulations, we see that they appear qualitatively similar. After a spin-up period, periodic events nucleate at a depth of approximately 10 km, and accumulate ∼3.5 m of slip at Earth’s surface. Comparing these to Fig. 2(d), where events nucleate further updip, and accumulate only ∼3 m slip with each event, we note that a decrease in either μ₁ or μ₃ (or both) increases the recurrence interval and thus the amount of slip during each earthquake. These results suggest that μ*, rather than absolute values of μ₁ and μ₃ determine model outcomes, including recurrence interval and nucleation depth.

Figure 2.

Cumulative slip plotted in blue every 5-a during the interseismic periods and in red every second during quasi-dynamic rupture for (a) an HTI simulation with μ₁ = 34 GPa, μ₃ = 16.95 GPa, (b) a VTI simulation with μ₁ = 16.95 GPa, μ₃ = 34 GPa, (c) an isotropic reference with μ* = 24 GPa and (d) an isotropic reference case with μ₁ = μ₃ = 34 GPa.

Since results remain similar for all cases with μ* = 24 GPa, we hypothesize that μ* predominantly determines emergent behaviour. To explore this, Fig. 3 presents several cases with μ* = 18 GPa and differing combinations of μ₁ and μ₃. Here we have an alternating sequence of large and small events nucleate at a depth of ∼12 km, a result seen by Lapusta & Rice (2003) for a decrease in h* (obtained by reducing the value of D_c rather than in shear modulus). This behaviour persists with stronger HTI anisotropy as seen in Fig. 3(a), with μ₁ = 36 GPa and μ₃ = 9 GPa. Interestingly, however, with the anisotropy reversed, the VTI simulation with μ₁ = 13.5 GPa and μ₃ = 24 GPa exhibits a more complex sequence of large, small, and medium sized events, as seen in Fig. 3(d). This VTI simulation compares most closely with the isotropic reference case corresponding to μ* = 18 GPa, as seen in Fig. 3(c). Comparing all results shown in Fig. 3 we observe that similar μ* values can generate quite different behaviours. We refer to these different types of behaviours as period two (in which large and small events emerge) and period three (where small, large, and medium sized events emerge).

Figure 3.

Cumulative slip plotted in blue every 5-a during the interseismic periods and in red every second during quasi-dynamic rupture for (a) an HTI simulation with μ₁ = 24 GPa, μ₃ = 13.5 GPa, (b) a second HTI simulation with μ₁ = 36 GPa, μ₃ = 9 GPa, (c) a VTI simulation with μ₁ = 13.5 GPa, μ₃ = 24 GPa and (d) an isotropic reference case with μ₁ = μ₃ = 18 GPa.

7.3 Aseismic Transients

In this parameter study we find that many simulations host transient events, where small increases in maximum slip rate emerge between large events. Aseismic transients are observed in other numerical simulations of earthquake sequences and are of interest to the broader community because they might indicate an impending large event (Lapusta & Liu 2009; Noda et al. 2013; Noda & Hori 2014).

In a study of transient events that emerge in earthquake cycle simulations within an isotropic medium, Noda & Hori (2014) find that A/B is a significant parameter that controls interseismic behaviour within a seismogenic patch, where A = aσ_n is the direct effect and B = bσ_n is the evolution effect in the rate-and-state friction law. They deduce that for A/B ≥ 0.6, aseismic transients events emerge when creep penetrates sufficiently far into the velocity-weakening zone to violate linear stability, a length scale given by h*, but before the creeping region can accommodate dynamic rupture, a length scale referred to as the nucleation size, which scales directly with the shear modulus μ. If we interpret this length scale for our anisotropic problems to scale with μ* we can contextualize our results in terms of these findings.

In our simulations the rate-and-state parameters a and b and the normal stress σ_n are fixed for all simulations. Thus A and B are fixed, with A/B = 0.75 within the seismogenic zone. We define an aseismic transient as an event where the maximum slip velocity (taken over the whole fault) climbs above a threshold of 10⁻⁹ m s⁻¹ but remains below the threshold of 1 mm s⁻¹ (which we define to be the threshold for coseismic speeds). In Fig. 4(a) we plot the time series of maximum slip velocity for the simulations from Fig. 2. The first three simulations from Fig. 2 share the same μ* value and have close to identical recurrence intervals of about 110 yr, but have some subtle differences. We observe that the HTI case (yellow) and the isotropic reference case (blue) both host aseismic transients, where small increases in maximum slip velocity occur approximately 20 yr before each large event. The VTI case with the same μ* value (red), however, does not host transients. In Fig. 4(b) we plot the nucleation zones for the simulations in Fig. 2, where fault shear stress is plotted as function of depth at the moment when rupture accelerates to coseismic speeds. The nucleation zone width we approximate numerically by computing the distance between the stress peaks surrounding the accelerating slip patch as in Rubin & Ampuero (2005). Smaller peaks further updip correspond to where interseismic creep has penetrated into the velocity weakening zone. The first three simulations from Fig. 2 all have the same μ*, and thus the same values of h* and nucleation size. And yet aseismic transients exist for only two of the simulations. For these simulations we observe that interseismic creep penetrates further updip for larger ratios R = μ₁/μ₃, before nucleation takes place, and numerical measurements show quite similar nucleation lengths of ∼2 km; see Fig. 4(b). The isotropic reference case in purple (with a larger μ*) nucleates higher updip and has a nucleation length of ∼3 km.

Figure 4.

(a) Maximum slip velocity time-series and (b) fault shear stress plotted against depth for simulations from Fig. 2 with μ* = 24 GPa. The isotropic case is plotted with a dashed line.

Fig. 5 shows the time series of maximum slip velocity, and the nucleation zones for all four simulations from Fig. 3. In Fig. 5(a), we observe that all four simulations host aseismic transients. The HTI simulations (blue and red) each host an aseismic transient before a large event, while the isotropic reference case with the same μ* (yellow) and the VTI case (purple) host transients before medium events. Examining the recurrence of events we see that for the HTI simulations, the interval before a large event is ∼161 yr, and ∼65 yr before a small event. The isotropic reference and VTI simulations, on the other hand, show some differences. The interval before a small event in the period three cycle, is ∼84 yr for both simulations. However, the interval before medium and large events differs for these two simulations: ∼106 and ∼178 yr, respectively, for the isotropic reference but only ∼96 and ∼184 yr (respectively) for the VTI case. We note that since the isotropic reference case hosts larger aseismic transients than the VTI, the difference in recurrence times between events may be due to this feature. Fig. 5(b) shows fault shear stress plotted against fault depth at the moment that maximum velocity reaches the threshold for coseismic slip. Residual stress from small sub-surface events lingers at ∼5 km depth for all simulations. In addition, between ∼5 and ∼6 km depth we see a spike in shear stress left over from the aseismic transient—we refer to in the figure as a failed rupture, or rupture that failed to reach coseismic speeds.

Figure 5.

(a) Maximum slip velocity time-series and (b) fault shear stress plotted against depth for simulations from Fig. 3 with μ* = 18 GPa. The isotropic case is plotted with a dashed line.

Just as in Fig. 4(b), larger ratios R = μ₁/μ₃ correspond to interseismic creep penetration further updip. Nucleation zone length and location is quite similar for all events, with length of ∼1.4 km, and location between ∼10 and ∼12 km depth. In this parameter regime, therefore, the emergence of aseismic transients cannot be attributed solely to A, B and μ*. The ratio R seems to determine how far creep penetrates updip before nucleation takes. The results we have presented thus far begin to reveal a more complex relationship between R, the nucleation zone, the presence of transients, and interseismic creep penetration. However, it is challenging to isolate the influence of any one of these factors.

7.4 Periodicity in Parameter Space

To better understand what determines period one, two or three behaviours, like those evidenced in Figs 2 and 3, we conducted a broad sweep of parameter values and list the results in Table 5. Descending the rows of Table 5 correspond to an increase in μ₁, while moving left to right corresponds to increasing values of μ₃. The value in the interior of each cell is the corresponding value of μ*. The colours of each cell correspond to the period of the simulation. Yellow indicates that the simulation has period one, red indicates period two, and blue indicates period three. Simulations with the same μ* value, for which differences in period occur when μ₁ and μ₃ are varied, are circled to highlight these differences. For example simulations with (μ₁, μ₃) = (30, 15) (row 12, column 4) and (15, 30) (row 4, column 12) share a μ* value of 21.21 but are period three and period one respectively. Similarly, μ* = 18 for simulations (μ₁, μ₃) = (24, 13.5) and (18, 18), but (μ₁, μ₃) = (24, 13.5) is period two while (μ₁, μ₃) = (18, 18) is period three. Bold cell value denote that the simulation hosts aseismic transients. White cells are simulations that are outside of the scope of this parameter study.

Table 5 reveals a bifurcation from period one to period two behaviour by decreasing μ*, with a complex boundary between regimes where period 3 behaviour emerges. The circles reveal that simulations with the same μ* can have quite different behaviours and periodicity, and the transition from italicized to bold fonts reveal that for most parameter values, a decrease in μ* corresponds to a transition from period one, to period one with aseismic transients, to period three with no aseismic transients, to period three with aseismic transients, to period two.

To ascertain more about the relationship between the ratio R = μ₁/μ₃, nucleation zone, and aseismic transients, we examine column 10 of Table 5, that is the column where μ₃ = 24. Note that to explore the role of R without having to vary both μ₁ and μ₃ at once (which further complicates analysis), it is necessary to instead vary μ*. In Fig. 6 we present the maximum slip velocity time series, and nucleation zones for a representative subset of the period one simulations from column 10, with a range of ratios of μ₁/μ₃. Based on the results presented thus far, we expect μ₁ to dictate how far updip interseismic creep is able to penetrate. As such, we would expect the extent of creep to be farthest updip for the highest ratio of R (shown in teal for R = 1.25) and lowest for the smallest ratio (shown in blue for R = 0.75). However, this is not entirely the case. Although the R = 0.75 case has the smallest nucleation zone and the smallest μ₁, interseismic creep is able to penetrate further up than the cases shown in red and yellow. Moreover, creep penetrates about as far updip for with larger μ₁/μ₃ ratios, namely, the red and yellow curves corresponding to R = 1 and R = ≈0.87, respectively. We attribute this to the influence of the interseismic transients present in both the R = 1 (yellow) and R = 0.75 (blue) simulations, (as evidenced in Fig. 6a). These results suggest that while R seems to govern nucleation zone size, and largely dictates how far updip interseismic creep can penetrate, the presence and magnitude of aseismic transients also plays a role.

Figure 6.

Maximum slip velocity and fault shear stress profiles shown for a representative sample of the simulations in Table 5 where μ₁ is allowed to vary and μ₃ is held constant at 24 GPa. The isotropic case is plotted with a dashed line.

To further explore aseismic transients events in other parts of parameter space, we looked more closely at some results from Table 5. In Fig. 7 we present plots of slip profiles for four simulations where μ₁ or μ₃ is set to 20.78 GPa and the other parameter takes on the values 11.972 GPa or 13.5 GPa. Figs 7(a) and (b) show a VTI and an HTI simulation corresponding to μ* = 15.77 GPa. For these parameter values, the same μ* can generate quite different results. Both simulations generate sequences of large and small events, but the large events in the VTI simulation have a longer recurrence interval, with more slip occurring with each event. Figs 7(c) and (d), however, show a VTI and an HTI simulation both with μ* = 16.75 GPa, where qualitatively similar sequences of events are generated. Figs 7(a) and (d) can be used to observe the effect of increasing μ₁, while keeping μ₃ fixed, while Figs 7(b) and (c) can be used to observe the effect of fixing μ₁ and varying μ₃.

Figure 7.

Cumulative slip plotted in blue every 5-a during the interseismic periods and in red every second during quasi-dynamic rupture for (a) a VTI simulation with μ₁ = 11.97 GPa, μ₃ = 20.78 GPa, (b) an HTI simulation with μ₁ = 20.78 GPa, μ₃ = 11.97 GPa, (c) an HTI simulation with μ₁ = 20.78 GPa, μ₃ = 13.5 GPa and (d) a VTI simulation with μ₁ = 13.5 GPa and μ₃ = 20.78 GPa.

In Fig. 8 the maximum slip velocity time series of all four simulations from Fig. 7 is plotted over a period of 350 yr. We see that all but one of the simulations hosts aseismic transients. The VTI simulation (in red) with μ₁ = 20.78 GPa, μ₃ = 11.97 GPa does not host aseismic transients, while its HTI counterpart (shown in yellow) does. We examine the role of shear stress in more detail in Fig. 8(b), where the nucleation zones are shown for all four simulations. Higher ratios of μ₁/μ₃ again correspond to further penetration updip of the interseismic creeping region, but does not correspond to whether or not transient events emerge, complicating the explanation given for transient presence and magnitude for the results in Fig. 7.

Figure 8.

(a) Maximum slip velocity and (b) fault shear stress plotted against depth shows nucleation zones for simulations from Fig. 7.

7.5 Full Anisotropy

As a final study, we report on some results when considering full anisotropy, i.e. where μ₂ ≠ 0. Fig. 9 shows results for both μ* = 24 GPa and 18 GPa, which can be compared to the orthotropic results with similar μ* in Figs 2 and 3. We found that these simulations required quite long spin-up periods, thus we plot cumulative slip profiles here relative to a background slip profile, so that figures may be more easily visualized. That is, we plot slip profiles relative to the total slip accumulated during the spin-up period.

Figure 9.

Cumulative slip plotted in blue every 5-a during the interseismic periods and in red every second during quasi-dynamic rupture for fully anisotropic simulations. In (a) and (b) μ* = 24, while in (c) and (d) μ* = 18 GPa. All simulations are plotted relative to a background slip profile.

Fig. 9(a) shows results for (μ₁, μ₂, μ₃) = (18.39, 7, 34) GPa, corresponding to a fast wave direction of about 20° from fault parallel (within the vertical plane) which we will refer to as ”near-VTI”. Fig. 9(b) is the reverse of this orientation (34, 7, 18.39) with a fast wave speed 20 degrees from fault perpendicular, which we will refer to as ”near-HTI”. Both cases provide additional evidence that model behaviours with μ* = 24 GPa are qualitatively similar, despite differences in μ₁, μ₂ and μ₃. However, this result does not persist for other values of μ*. Fig. 9(c) shows results for a near-VTI simulation with μ* = 18 GPa, with (μ₁, μ₂, μ₃) = (14.56, 5, 24), and a fast-wave direction again about 20° from fault-parallel. Results are qualitatively similar to the VTI simulation shown in Fig. 3(d), where small, medium, and large events emerge. Fig. 9(d) shows results for an orientation reverse with (μ₁, μ₂, μ₃) = (24, 5, 14.56), with results more similar to the HTI simulation shown in Fig. 3(c) where only large and small events emerge.

In Fig. 10 we examine the maximum slip velocity time series and nucleation zone profiles for simulations in Fig. 9. The near-HTI simulations (red and purple curves) as well as the one with near-VTI (yellow) all host transients, whereas the near-VTI simulation in blue does not. The near-VTI simulation in yellow hosts three transients: one in the recurrence interval leading up to a small event and two more in the recurrence interval preceding a large event. In the latter interval a smaller transient occurs ∼38 yr before a large event, and then another much larger in magnitude transient occurs ∼8 yr before a large event. In Fig. 10(b) we see that the large transient corresponds to a failed event that was able to partially propagate up and down the fault, destabilizing the patch of the fault between ∼3.3 and 7.3 km depth, but was unable to nucleate to coseismic speeds. We suspect the large event corresponding to the near-VTI simulation (in yellow) successfully nucleates further updip than the large event in the near-HTI simulation (in purple, a simulation with the same μ*) as a result of the residual stress from the large aseismic transient ∼8 yr prior. These results complicate our findings in the orthotropic parameter study, where we observed that transients appear to allow creep to penetrate further updip. For this fully anisotropic parameter regime, it appears as though transients allow nucleation to occur in the presence of less updip creep. These seemingly conflicting results may be reconciled if we allow for the possibility that the timing within the interseismic period of an aseismic transient also plays a role. As such, the temporal location of aseismic transients should be investigated further in future work.

Figure 10.

(a) Maximum slip velocity and (b) fault shear stress plotted against depth shows nucleation zones for simulations from Fig. 9 with full anisotropy.

7.6 Surface Velocity Profiles

In this section we report on surface velocity profiles which could be linked to observables measured at Earth’s surface across a fault. Fig. 11 shows plots of surface velocity as a function of off-fault distance (y) for the results shown so far. Panels (a)–(c) correspond to orthotropic scenarios and panel (d) [and corresponding zoom in panel (e)] shows surface velocities for the fully anisotropic scenarios. Dashed and solid contours correspond respectively to 25 and 95 per cent through the recurrence interval preceding a large event for each simulation. Fig. 11(a) shows the three similar sequences corresponding to the same μ* value (blue, red, yellow). In both time instances, the recurrence interval we observe an increase in strain (∂u/∂y) corresponding to a decrease in μ₁, which is to be expected if shear stress is to remain constant. Strain increases with decreasing μ₁ also for the cases shown in panels (b) and (c). For the fully anisotropic simulations (d) and zoom (e), this feature persists when comparing those simulations with equivalent μ*. We include the zoom (e) to illustrate the full anisotropy can allow for small amounts of surface creep during the interseismic period, due to the fact that the ∂u/∂z component of strain is non-negligible and contributes to the fault shear stress.

Figure 11.

Surface velocity profiles plotted before a large event, after each simulation is out of its spin-up cycle, with dashed lines 25 per cent of the way into the recurrence interval and solid lines 95 per cent of the way into the recurrence interval. Panels (a)–(c) correspond to the orthotropic simulations from Figs 2 to 7, while panel (d) and corresponding zoom (e) correspond to the fully anisotropic simulations shown in Fig. 9. The Isotropic cases are indicated with a small I placed in the legend next to μ₁.

8 DISCUSSION AND FUTURE WORK

We note that in developing the method for the 2-D antiplane problem, we have inherently limited the possible directions for wave-propagation. In fully anisotropic simulations, like the ones in Section 7.5, waves propagate at oblique angles to the fault which cannot be reliably observed in practice. However, real-world observations of fault parallel fast waves, like those in the VTI simulations above, have been made (Stuart et al. 2002).

Although the antiplane setting has been explored exclusively in this work, extensions to the in-plane setting will also be insightful due to the richer range of wave speeds. The computational framework for earthquake cycles in 2-D plane strain has been developed for the isotropic setting in Erickson and Day (2016) and is readily extendable to anisotropic material. An anisotropic response will generate perturbations in normal stress, even for monomaterial interfaces, due to additional components of strain that are incorporated and may lead to the emergence of Weertman pulses seen on bimaterial interfaces (Weertman 1980; Andrews & Ben-Zion 1997; Anooshehpoor & Brune 1999)

The results of our parameter studies, show that the size and location of the nucleation zone for a simulation, is influenced not just by the ratio μ₁/μ₃, but also by the presence and magnitude of aseismic transients. This suggests that both play a role in how far updip interseismic creep may penetrate. However, the emergence of more complicated aseismic transients in the fully anisotropic simulations leaves questions to be explored about the relevance of the temporal location of these transients. Additionally, relationships and interactions between the residual stresses from subsurface events, failed ruptures near these and aseismic transients should be explored with further studies.

9 CONCLUSION

We have extended the computational framework developed in Erickson & Dunham (2014) and adapted it to study earthquake cycles in anisotropic media. The off-fault volume is discretized with finite difference operators satisfying a summation-by-parts rule, with weak enforcement of boundary conditions, which leads to a provably stable formulation. Rate-and-state friction is enforced along the fault, and sequences of earthquakes are generated by displacing the remote boundary at a slow plate rate. We tested our numerical scheme by applying it to a suitable manufactured solution and ensuring it achieved the expected order of convergence.

In the parameter studies in this work, we found that anisotropy influences the recurrence interval, periodicity, emergence of transients, nucleation zone size and depth, and extent of interseismic creep penetration. We found that choices for μ₁/μ₃ can cause simulations with the same μ* value to exhibit quite different behaviour, and uncovered a complex boundary between the period one and period two solutions that naturally arise as one decreases h*. We found that this boundary often exhibits period three behaviour and gives rise to simulations that host aseismic transients. We additionally found that in period two and three solutions, residual stresses from subsurface events appear, and for some simulations failed rupture often occurs near these residual stresses ahead of a coseismic event.

ACKNOWLEDGEMENTS

MBM and BAE were supported through the NSF under Award No. EAR-1547603 and by the Southern California Earthquake Center. SCEC is funded by NSF Cooperative Agreement EAR-0529922 and USGS Cooperative Agreement 07HQAG0008 (SCEC contribution number 9205). JEK was supported under NSF Award No. EAR-1547596. Thanks to Elizabeth Cochran for useful discussion. We would also like to acknowledge the two reviewers and the editor Eiichi Fukuyama whose comments and suggestions much improved this work. We would also like to thank Alice-Agnes Gabriel, the anonymous reviewer, and editor Eiichi Fukuyama whose comments and suggestions much improved this work.

APPENDIX A: RESOLUTION OF THE COHESIVE ZONE

We show our simulations are well-resolved for a period three isotropic problem with μ = 18 GPa. Fig. A1 shows cumulative slip profiles, where L_b is resolved with over 5 grid points on the left and over 10 on the right. We see differences only during the spin-up cycle of each simulation, after which both settle into period three behaviour that is qualitatively similar.

Figure A1.

Cumulative slip plotted in blue every 5-a during the interseismic periods and in red every second during quasi-dynamic rupture for two, period three, isotropic simulations, where all parameters are held constant. In panel (a) the number of grid points is set to 1165 to ensure that cohesive zone L_b is resolved with over 5 grid points, while in panel (b) the grid points are more than doubled to 2500 and L_b is resolved with over 10 grid points.

APPENDIX B: EDGE EFFECTS

Truncating the off-fault computational domain, L_y, at 72 km causes large events in some of the multiperiod simulations in our parameter study to nucleate at an artificially shallow depth of around 5 km. We suspect that this is due to the influence of finite distance to the remote boundary where loading is enforced. We doubled the domain size for one such simulation, a period two simulation with μ₁ = 36 GPa, and μ₃ = 9 GPa. Increasing the computational domain size leads to large events that nucleate farther downdip, closer to 12 km depth. In Fig. B1, we compare the cumulative slip plots with a domain size of L_y = 72 km on the left, and L_y = 144 km on the right. It is worth noting that this edge effect, when it occurs, appears in HTI orthotropic simulations where μ₁ > μ₃, but not in the VTI counterpart (with the same μ*) for which the values of μ₁ and μ₃ are reversed. Simulations where these edge effect occurred were re-run with a doubled computational domain, and replaced throughout the paper.

Figure B1.

Cumulative slip plotted in blue every 5-a during the interseismic periods and in red every second during quasi-dynamic rupture for two orthotropic simulations with μ₁ = 36 and μ₂ = 9. In panel (a) L_y = 72 km. In panel (b) L_y = 144 km, that is, the off-fault computational domain is doubled.