Introduction

The earthquake cycle involves many processes having different time and length scales, such as interseismic stress accumulation, nucleation and propagation of dynamic rupture, and postseismic deformation. In modeling each of these serial processes, the initial condition is essential for determining its outcome, which in turn gives the initial condition for the modeling of the following process. Therefore, to investigate the behavior of a fault that accommodates repeating earthquakes, these processes must be simulated in a consistent manner within a single framework, without changing the governing equations artificially. Dynamic earthquake sequence simulation (e.g., Rice and Ben-Zion 1996; Lapusta et al. 2000) enables us to numerically simulate the behavior of a fault embedded in a continuum such that both slow interseismic processes and rapid coseismic processes, in which the inertial, wave-mediated effect matters, are taken into account. Dynamic earthquake sequence simulation has been used to investigate various important problems in earthquake physics, such as the effects of depth-dependent frictional properties (e.g., Lapusta et al. 2000; Shimamoto and Noda 2014), brittleness of the seismogenic part of a fault (e.g., Lapusta and Rice 2003; Chen and Lapusta 2009; Noda et al. 2013; Barbot 2019), drastic weakening of a fault at coseismic slip rate (e.g., Noda and Lapusta 2013), and off-fault viscous deformation (e.g., Miyake and Noda 2019).

Because dynamic earthquake sequence simulation is computationally intense, a spectral boundary integral equation method (SBIEM) is often used for the simulation, even though the existing SBIEMs are applicable only to problems with simple geometry. In applying a boundary integral equation method, only the boundaries where the traction needs to be evaluated and the sources of the stress change are distributed must be discretized. Thus, simulations having a greater on-fault resolution can be conducted with restricted computational resources, compared with methods based on discretization of the volume. In addition, the fast Fourier transform (FFT) algorithm can be used to reduce the computational cost of the spectral method.

Figure 1a schematically shows a typical problem. A simulated part of a fault governed by a rate- and state-dependent friction law (Dieterich 1979; Ruina 1983) is loaded by background uniform slip rate \(V_{{{\text{pl}}}}\). It has length \(L\), typically has rate-strengthening and rate-weakening parts, and the slip difference from the uniform background level \(\delta ^{\prime}\) distributes within it. Earthquake ruptures occur in the rate-weakening patches. Although the spectral method is very efficient, Fourier series representation of a Green’s function in an infinite elastic space cause periodic replication if the wavenumber domain is discretized by using regularly spaced grid points (Fig. 1b). To obtain a good approximation of fault behavior in an infinite medium, the distance \(\lambda\) between the periodic boundaries should be much larger than the length \(L\) of the system of interest. For example, Lapusta et al. (2000) used \(\lambda = 4L\).

Fig. 1
figure 1

a A schematic diagram of a typical problem solved with dynamic earthquake sequence simulation. A region with source of traction change (slip difference from a reference \(\delta ^{\prime}\)) distributes in a region of length \(L\), which is governed by a rate- and state-dependent friction law. b Spectral representation and discretization with a regular mesh in the wavenumber domain causes artificial periodic boundaries every \(\lambda\). Computational domains required in the present method and in the conventional method are indicated by dashed rectangles in a and b, respectively

Full size image

Cochard and Rice (1997) developed a novel method to remove the periodic boundary condition from an SBIEM for simulations of dynamic rupture. For slip δ distributed over a fault of length L, they prepared an additional zero-padded region of length L outside the source region. Thus, they expressed the distribution of the source as a Fourier series for a domain of length 2L, and derived an expression for the integration kernel in the wavenumber domain that gave the stress change without periodic replication, provided that the collocation point was within the source region. In the present study, a method for dynamic earthquake sequence simulation using a new SBIEM that partly incorporates their method is proposed.

It turns out that calculation of the dynamic part of the integration kernel for a dynamic earthquake sequence simulation incorporating the formulation by Cochard and Rice (1997) using a velocity formulation (e.g., Rice and Ben-Zion 1996; Lapusta et al. 2000) requires evaluation of a two-dimensional numerical integral. Because this calculation takes an unnecessarily long time, the method by Cochard and Rice (1997) is applied only to the static part of the integration kernel, which has an infinitely long-range effect with a spatially compact source distribution.

Methodology

SBIEM for dynamic earthquake sequence simulation

A two-dimensional problem of a planar fault embedded in a linearly elastic infinite body is considered. The traction on the fault can be written as follows (Cochard and Madariaga 1994; Geubelle and Rice 1995):

$$\begin{array}{*{20}c} {\tau \left( {x,t} \right) = \tau _{0} \left( {x,t} \right) + \phi - \eta V^{\prime}\left( {x,t} \right)~,} \\ \end{array}$$
(1)

where \(x\) is the coordinate along the fault, \(t\) is time, \(\tau\) is the traction on the fault, \(\tau _{0}\) is the traction realized in a reference state (e.g., no relative motion, steady-state solution, or uniform slip rate), \(\phi\) is the wave-mediated stress transfer due to previous fault motion relative to the reference, and \(\eta\) is a proportionality constant between instantaneous changes in \(\tau\) and the rate of fault motion relative to the reference state \(V^{\prime}\). Calculation of \(\phi\) requires spatio-temporal convolution. The last term on the right-hand side represents the impedance effect, which is sometimes called radiation damping, and \(\eta\) depends on the mode of the rupture,

$$\begin{array}{*{20}c} {\eta = \left\{ {\begin{array}{*{20}c} {\rho c_{\text{s}} /2} & {\left( {{\text{Modes II and III}}} \right)} \\ {\rho c_{\text{p}} /2} & {\left( {{\text{Mode I}}} \right)} \\ \end{array} } \right.,} \\ \end{array}$$
(2)

where \(\rho\) is the density and \(c_{{\text{s}}}\) and \(c_{{\text{p}}}\) are the S-wave speed and P-wave speed, respectively.

If the fault motion is suddenly stopped at time \(t\), then the elastic waves propagate away from the fault and the static solution of the 2-dimensional problem is realized asymptotically. The traction change in this static solution is denoted as \(\phi _{{{\text{st}}}} \left( t \right)\), which can be expressed by the spatial convolution of a static Green’s function and the distribution of the relative displacement difference from the reference state at time \(t\): \(\delta ^{\prime}\left( {x,t} \right)\). The dynamic part of \(\phi\), \(\phi _{{{\text{dyn}}}}\), is defined as follows (Rice and Ben-Zion 1996):

$$\begin{array}{*{20}c} {\phi _{{{\text{dyn}}}} = \phi - \phi _{{{\text{st}}}} .} \\ \end{array}$$
(3)

In the displacement formulation, \(\phi\), is expressed as

$$\begin{array}{*{20}c} {\phi \left( {x,t} \right) = \mathop \smallint \limits_{{ - \infty }}^{t} {\text{d}}t^{\prime}\mathop \smallint \limits_{{ - \infty }}^{\infty } {\text{d}}x^{\prime}\left[ {G\left( {x - x^{\prime},t - t^{\prime}} \right)\delta ^{\prime}\left( {x^{\prime},t^{\prime}} \right)} \right]~,} \\ \end{array}$$
(4)

where \(G\) is a Green’s function representing the traction change due to spatio-temporally localized relative displacement. Fourier transformation in terms of \(x\) leads to

$$\begin{array}{*{20}c} {\mathop \smallint \limits_{{ - \infty }}^{\infty } \phi \left( {x,t} \right)e^{{ - \text{i}kx}} {\text{d}}x = \Phi \left( {k,t} \right) = \mathop \smallint \limits_{{ - \infty }}^{t} K\left( {t - t^{\prime}} \right)D\left( {k,t^{\prime}} \right){\text{d}}t^{\prime}} \\ { = - \frac{{\mu c_{\text{s}} k^{2} }}{2}\mathop \smallint \limits_{{ - \infty }}^{t} C\left( {\left| k \right|c_{\text{s}} \left( {t - t^{\prime}} \right)} \right)D\left( {k,t^{\prime}} \right){\text{d}}t^{\prime},} \\ \end{array}$$
(5)

where \({\text{i}} = \sqrt { - 1}\), \(k\) is the angular wavenumber and \(\mu\) is the shear modulus. \(K\) is an integration kernel in the angular wavenumber domain in the displacement formulation, and \(D\) is the Fourier transform of \(\delta ^{\prime}\):

$$\begin{array}{*{20}c} {D\left( {k,t} \right) = \mathop \smallint \limits_{{ - \infty }}^{\infty } \delta ^{\prime}\left( {x,t} \right)e^{{ - \text{i}kx}} {\text{d}}x.} \\ \end{array}$$
(6)

\(C\) is a function appearing in the kernel and depends on the mode of the rupture, as described in Appendix (Geubelle and Rice 1995). It was calculated by polynomial approximation (Morrissey and Geubelle 1997) and integrated numerically in previous studies (e.g., Lapusta et al. 2000). In the simulation code used in the present paper, it was calculated using special functions (see Appendix) implemented in the Python routine scipy.special. To extract the static term, integration by parts of Eq. (5) yields

$$\begin{array}{*{20}c} {\Phi \left( {k,t} \right) = - \frac{{\mu c_{\text{s}} k^{2} }}{2}\left[ {\mathop \smallint \limits_{{ - \infty }}^{t} C\left( {\left| k \right|c_{\text{s}} \left( {t - t^{\prime}} \right)} \right){\text{d}}t^{\prime}} \right]D\left( {k,t} \right)} \\ { + \frac{{\mu c_{\text{s}} k^{2} }}{2}\mathop \smallint \limits_{{ - \infty }}^{t} \left[ {\mathop \smallint \limits_{{ - \infty }}^{{t^{\prime}}} C\left( {\left| k \right|c_{\text{s}} \left( {t - t^{\prime\prime}} \right)} \right){\text{d}}t^{\prime\prime}} \right]\dot{D}\left( {k,t^{\prime}} \right){\text{d}}t^{\prime}~,} \\ \end{array}$$
(7)

where \(\dot{D}\) is the Fourier transform of \(V^{\prime}\). It has been assumed that \(D\left( {k,t} \right)\) is zero at sufficiently small \(t\). The first term on the right-hand side of Eq. (7) represents the Fourier transform of \(\phi _{{{\text{st}}}}\),

$$\begin{array}{*{20}c} {\mathop \smallint \limits_{{ - \infty }}^{\infty } \phi _{{{\text{st}}}} \left( {x,t} \right)e^{{ - \text{i}kx}} {\text{d}}x = \Phi _{{{\text{st}}}} \left( {k,t} \right) = - \frac{{\mu c_{\text{s}} k^{2} }}{2}\left[ {\mathop \smallint \limits_{{ - \infty }}^{t} C\left( {\left| k \right|c_{\text{s}} \left( {t - t^{\prime}} \right)} \right){\text{d}}t^{\prime}} \right]D\left( {k,t} \right).} \\ \end{array}$$
(8)

After executing the integral, Eq. (8) can be written as

$$\begin{array}{*{20}c} {\Phi _{{{\text{st}}}} \left( {k,t} \right) = K_{{{\text{st}}}} \left( k \right)D\left( {k,t} \right),} \\ \end{array}$$
(9)

where \(K_{{{\text{st}}}}\) is the static kernel in the angular wavenumber domain,

$$\begin{array}{*{20}c} {K_{{{\text{st}}}} \left( k \right) = - \frac{{\mu ^{\prime}\left| k \right|}}{2}~.} \\ \end{array}$$
(10)

\(\mu ^{\prime}\) depends on the direction of the relative displacement: \(\mu ^{\prime} = \mu\) in mode III problems, and \(\mu ^{\prime} = \mu /\left( {1 - \nu } \right)\) in mode I and mode II problems, where \(\nu\) is the Poisson’s ratio. The second term on the right-hand side in Eq. (7) represents the dynamic contribution,

$$\begin{array}{*{20}c} {\mathop \smallint \limits_{{ - \infty }}^{\infty } \phi _{{{\text{dyn}}}} \left( {x,t} \right)e^{{ - \text{i}kx}} \text{d}x = \Phi _{{{\text{dyn}}}} \left( {k,t} \right) = \mathop \smallint \limits_{{ - \infty }}^{t} K_{{{\text{dyn}}}} \left( {k,t - t^{\prime}} \right)\dot{D}\left( {k,t^{\prime}} \right){\text{d}}t^{\prime},} \\ \end{array}$$
(11)

where the dynamic kernel is

$$\begin{array}{*{20}c} {K_{{{\text{dyn}}}} \left( {k,t} \right) = \frac{{\mu \left| k \right|}}{2}\mathop \smallint \limits_{{\left| k \right|c_{\text{s}} t}}^{\infty } C\left( s \right){\text{d}}s.} \\ \end{array}$$
(12)

\({{\Phi }}\) can therefore be expressed as

$$\begin{array}{*{20}c} {\Phi \left( {k,t} \right) = K_{{{\text{st}}}} \left( k \right)D\left( {k,t} \right) + \mathop \smallint \limits_{{ - \infty }}^{t} K_{{{\text{dyn}}}} \left( {k,t - t^{\prime}} \right)\dot{D}\left( {k,t^{\prime}} \right){\text{d}}t^{\prime}.} \\ \end{array}$$
(13)

As discussed by Rice and Ben-Zion (1996) and demonstrated by Lapusta and Liu (2009), \(K_{{{\text{dyn}}}}\) decays to zero at large \(t - t^{\prime}\) so that the effect of an old source becomes negligible. Therefore, we can truncate the time integration in Eq. (13) by assuming a time window \(t_{{\text{w}}}\) within which the inertial effect is fully accounted for:

$$\begin{array}{*{20}c} {\Phi \left( {k,t} \right) \approx \Phi _{\text{w}} \left( {k,t;t_{\text{w}} } \right) = K_{{{\text{st}}}} \left( k \right)D\left( {k,t} \right) + \mathop \smallint \limits_{{t - t_{\text{w}} }}^{t} K_{{{\text{dyn}}}} \left( {k,t - t^{\prime}} \right)\dot{D}\left( {k,t^{\prime}} \right){\text{d}}t^{\prime}.} \\ \end{array}$$
(14)

This calculation only requires a history with a finite length and thus enables us to simulate a sequence of earthquakes with a finite amount of computer memory.

Removal of periodic boundaries by the method of Cochard and Rice (1997)

In the angular wavenumber domain, the spatial convolution becomes multiplication for each \(k\). The use of a FFT reduces the computational cost and makes the SBIEM one of the most efficient methods for modeling of the fault behavior. However, discretization of the wavenumber domain with a regular mesh with interval \(\Delta k\) results in periodic boundaries at every \(\lambda = 2\pi /\Delta k\) (Fig. 1b). Cochard and Rice (1997) proposed a novel trick for removing the periodic boundary condition: for a function \(f\) having non-zero values only within \(\left[ { - L/2,L/2} \right]\),

$$\begin{array}{*{20}c} {\mathop \smallint \limits_{{ - \infty }}^{\infty } f\left( x \right){\text{d}}x = \mathop \smallint \limits_{{X - L}}^{{X + L}} f\left( x \right){\text{d}}x~~~{\text{for}}~X \in \left[ { - L/2,L/2} \right]~.} \\ \end{array}$$
(15)

Because the discretization of \(k\) introduces a periodic boundary, Cochard and Rice (1997) kept the domain of \(k\) continuous, but they used a different spectral expression for \(\delta ^{\prime}\). For \(\delta ^{\prime}\) with a non-zero value within a region of size \(L\), \(\left[ { - L/2,L/2} \right]\), they zero-padded outside it and expressed \(\delta ^{\prime}\) as a Fourier series for a function defined in \(\left[ { - L,L} \right]\),

$$\begin{array}{*{20}c} {\delta ^{\prime}\left( {x,t} \right) = \mathop \sum \limits_{{n = - \infty }}^{\infty } D_{n} \left( t \right)e^{{2\pi \text{i}nx/2L}} ,} \\ \end{array}$$
(16)

where \(n\) is the integer for the number of waves in \(\left[ { - L,L} \right]\). By plugging Eq. (16) into Eq. (6) and using Eq. (15), \(D\) can be expressed as

$$\begin{array}{*{20}c} {D\left( {k,t} \right) = \left( {\mathop \sum \limits_{{n = - \infty }}^{\infty } D_{n} \left( t \right)\frac{{2\sin \left( {kL - n\pi } \right)}}{{k - n\pi /L}}e^{{2\pi \text{i}nX/2L}} } \right)e^{{ - \text{i}kX}} } \\ \end{array}$$
(17)

for \(X \in \left[ { - L/2,L/2} \right]\). Inverse Fourier transformation of Eq. (5) with respect to \(k\) with Eq. (17) and evaluation of \(\phi\) at \(x = X\) results in

$$\begin{array}{*{20}c} {\phi \left( {x,t} \right) = \mathop \sum \limits_{{n = - \infty }}^{\infty } \Phi _{n} \left( t \right)e^{{2\pi \text{i}nx/2L}} ,} \\ \end{array}$$
(18)

where

$$\begin{array}{*{20}c} {\Phi _{n} \left( t \right) = - \frac{{\mu c_{\text{s}} }}{2}\left( {\frac{{\pi n}}{L}} \right)^{2} \mathop \smallint \limits_{{ - \infty }}^{t} C_{n}^{'} \left( {t - t^{\prime}} \right)D_{n} \left( {t^{\prime}} \right){\text{d}}t^{\prime}~,} \\ \end{array}$$
(19)
$$\begin{array}{*{20}c} {C_{n}^{'} \left( t \right) = \frac{{\left( { - 1} \right)^{n} }}{\pi }\mathop \smallint \limits_{{ - \infty }}^{\infty } \frac{{\sin \left( {kL} \right)}}{{k - n\pi /L}}C\left( {\left| k \right|c_{\text{s}} t} \right){\text{d}}k.} \\ \end{array}$$
(20)

Note that this expression is correct only for \(x \in \left[ { - L/2,L/2} \right]\). As pointed out by Cochard and Rice (1997), the kernel is identical to that with periodic boundaries when the delay time is small, that is, before the elastic waves from the periodic replication reach the region where the traction is evaluated.

To demonstrate their method, Cochard and Rice (1997) applied it to an antiplane problem of a crack under instantaneous loading. In the present study, this technique is applied to remove the periodic boundary condition in a dynamic earthquake sequence simulation for all modes.

Application to a dynamic earthquake sequence simulation

To apply the technique of Cochard and Rice (1997) to earthquake sequence simulation (e.g., Lapusta et al. 2000), the displacement formulation (Eq. (19)) is converted to a velocity formulation by conducting integration by parts:

$$\begin{array}{*{20}c} {\Phi _{n} \left( t \right) = K_{{{\text{st}}}}^{*} \left( n \right)D_{n} \left( t \right) + \mathop \smallint \limits_{{ - \infty }}^{t} K_{{{\text{dyn}}}}^{*} \left( {n,t - t^{\prime}} \right)\dot{D}_{n} \left( {t^{\prime}} \right){\text{d}}t^{\prime},} \\ \end{array}$$
(21)

where the superscript * indicates that the static and dynamic kernels are derived by using the formulation of Cochard and Rice (1997) and expressed in terms of \(C^{\prime}\) as

$$\begin{array}{*{20}c} {K_{{{\text{st}}}}^{*} \left( n \right) = - \frac{{\mu c_{\text{s}} }}{2}\left( {\frac{{\pi n}}{L}} \right)^{2} \mathop \smallint \limits_{0}^{\infty } C_{n}^{'} \left( s \right){\text{d}}s} \\ \end{array}$$
(22)

and

$$\begin{array}{*{20}c} {K_{{{\text{dyn}}}}^{*} \left( {n,t} \right) = \frac{{\mu c_{\text{s}} }}{2}\left( {\frac{{\pi n}}{L}} \right)^{2} \mathop \smallint \limits_{t}^{\infty } C_{n}^{'} \left( s \right){\text{d}}s,} \\ \end{array}$$
(23)

respectively.

For mode III (antiplane) problems,

$$\begin{array}{*{20}c} {C\left( T \right) = C_{{{\text{III}}}} \left( T \right) = \frac{{\text{J}_{1} \left( T \right)}}{T}.} \\ \end{array}$$
(24)

In this case, the integral in Eqs. (20) and (22) can be executed analytically as,

$$\begin{array}{*{20}c} {K_{{{\text{st}}}}^{{{\text{III}}*}} \left( n \right) = - \frac{{n\mu }}{L}{\text{Si}}\left( {n\pi } \right)} \\ \end{array} ,$$
(25)

where \({\text{Si}}\) is the sine-integral function

$$\begin{array}{*{20}c} {{\text{Si}}\left( x \right) = \mathop \smallint \limits_{0}^{x} \frac{{\sin \left( t \right)}}{t}{\text{d}}t.} \\ \end{array}$$
(26)

For the other rupture modes, we have only to modify the shear modulus depending on the Poisson’s ratio. Therefore, the static kernels for all modes are compactly expressed as

$$\begin{array}{*{20}c} {K_{{{\text{st}}}}^{*} \left( n \right) = - \frac{{n\mu ^{\prime}}}{L}{\text{Si}}\left( {n\pi } \right)} \\ \end{array} ,$$
(27)

similarly to the case with a periodic boundary (Eq. (10)). Because \({\text{Si}}\left( \infty \right) = \pi /2\), \(K_{{{\text{st}}}}^{*} \left( n \right)\) approaches \(K_{{{\text{st}}}} \left( {2\pi n/2L} \right)\) when the wavenumber is large.

Expressions for the dynamic kernel involve two-dimensional integration (Eqs. (20) and (23)), which the author could not fully solve analytically. For mode III, \(K_{{{\text{dyn}}}}^{*} \left( {n,t} \right)\) can be written using an integral expression for \(C_{{{\text{III}}}}\) (Cochard and Rice 1997) as

$$\begin{array}{*{20}c} {K_{{{\text{dyn}}}}^{{{\text{III}}*}} \left( {n,t} \right) = \frac{{\mu c_{\text{s}} }}{\pi }\left( {\frac{{\pi n}}{L}} \right)^{2} \mathop \smallint \limits_{t}^{\infty } \mathop \smallint \limits_{0}^{{\psi _{1} }} \cos ^{2} \left( \psi \right)\cos \left( {\frac{{n\pi c_{\text{s}} s}}{L}\sin \left( \psi \right)} \right){\text{d}}\psi {\text{d}}s,} \\ \end{array}$$
(28)

where

$${\psi _1} = \left\{ {\begin{array}{*{20}{c}} {\pi /2}&{{c_{\text{s}}}t/L \le 1}\\ {\arcsin (L/{c_{\text{s}}}t)}&{{c_{\text{s}}}t/L > 1} \end{array}} \right..$$
(29)

However, numerical integration of Eq. (28) with the Python routine scipy.integrate.dblquad turned out to take an excessively long time. Cochard and Rice (1997) proposed an advantageous way to evaluate the integral with respect to \(\psi _{1}\) for small \(c_{s} t/L\left( { > 1} \right)\), by subtracting integral from \(\psi _{1}\) to \(\pi /2\) from integral from \(0\) to \(\pi /2\). There may be an efficient way to calculate Eq. (28) and its counterparts in different modes, but here I propose an alternative method to remove the periodic boundary condition without using Eq. (23). In the present method, the dynamic kernel can be calculated without numerical integration if special functions in scipy.special are used (Appendix).

The fundamental idea of earthquake sequence simulation is truncation of the temporal convolution by setting a time window \(t_{{\text{w}}}\) (Eq. (14)),

$$\begin{array}{*{20}c} {{{\Phi }}_{{{\text{w}}n}} \left( {t;t_{{\text{w}}} } \right) = K_{{{\text{st}}}}^{*} \left( n \right)D_{n} \left( t \right) + \mathop \smallint \limits_{{t - t_{{\text{w}}} }}^{t} K_{{{\text{dyn}}}}^{*} \left( {n,t - t^{\prime}} \right)\dot{D}_{n} \left( {t^{\prime}} \right){\text{d}}t^{\prime}.} \\ \end{array}$$
(30)

\(\phi\), after truncation of the temporal convolution, can be written in two ways:

$$\begin{array}{l} \phi \left( {x,t} \right) \approx {\phi _{\text{w}}}\left( {x,t;{t_{\text{w}}}} \right) = \frac{1}{{2\pi }}\int\limits_{ - \infty }^\infty {{\Phi _{\text{w}}}} \left( {k,t;{t_{\text{w}}}} \right){e^{\text{i}kx}}{\rm{d}}k\\ \;\quad \quad \; = \sum\limits_{n = - \infty }^\infty {{\Phi _{\text{w}n}}} \left( {t;{t_{\text{w}}}} \right){e^{2\pi \text{i}nx/2L}}. \end{array}$$
(31)

Let \(\phi _{{\text{w}}}\) be split into two parts.

$$\begin{array}{*{20}c} {\phi _{{\text{w}}} \left( {x,t;t_{{\text{w}}} } \right) = \phi _{{{\text{past}}}} \left( {x,t;t_{{\text{w}}} } \right) + \Delta \phi \left( {x,t;t_{{\text{w}}} } \right),} \\ \end{array}$$
(32)

where

$$\begin{array}{*{20}c} {\phi _{{{\text{past}}}} \left( {x,t;t_{{\text{w}}} } \right) = \frac{1}{{2{{\pi }}}}\mathop \smallint \limits_{{ - \infty }}^{\infty } K_{{{\text{st}}}} \left( k \right)D\left( {k,t - t_{{\text{w}}} } \right){\text{e}}^{{{\text{i}}kx}} {\text{d}}k} \\ { = \mathop \sum \limits_{{n = - \infty }}^{\infty } K_{{{\text{st}}}}^{*} \left( n \right)D_{n} \left( {t - t_{{\text{w}}} } \right){\text{e}}^{{2{{\pi \text{i}}}nx/2L}} ~} \\ \end{array}$$
(33)

is the static effect in \(\phi _{{\text{w}}}\) due to sources before \(t - t_{{\text{w}}}\). Because the time derivative of \(D\) is \(\dot{D}\), the \(\Delta \phi\) can be written as

$$\begin{array}{*{20}c} {\Delta \phi \left( {x,t;t_{{\text{w}}} } \right) = \frac{1}{{2{{\pi }}}}\mathop \smallint \limits_{{ - \infty }}^{\infty } {\text{d}}k{\text{e}}^{{{\text{i}}kx}} \mathop \smallint \limits_{{t - t_{{\text{w}}} }}^{t} {\text{d}}t^{\prime}\left[ {\left( {K_{{{\text{dyn}}}} \left( {k,t - t^{\prime}} \right) + K_{{{\text{st}}}} \left( k \right)} \right)\dot{D}\left( {k,t^{\prime}} \right)} \right]~,} \\ \end{array}$$
(34)
$$\begin{array}{*{20}c} {\Delta \phi \left( {x,t;t_{{\text{w}}} } \right) = \mathop \sum \limits_{{n = - \infty }}^{\infty } {\text{e}}^{{2{{\pi \text{i}}}nx/2L}} \mathop \smallint \limits_{{t - t_{{\text{w}}} }}^{t} {\text{d}}t^{\prime}\left[ {\left( {K_{{{\text{dyn}}}}^{*} \left( {n,t - t^{\prime}} \right) + K_{{{\text{st}}}}^{*} \left( n \right)} \right)\dot{D}_{n} \left( {t^{\prime}} \right)} \right]~.} \\ \end{array}$$
(35)

Because the static stress change extends infinitely far from the source, it is important to use the technique by Cochard and Rice (1997) for its calculation. Therefore, for \(\phi _{{{\text{past}}}}\) the expression using \(K_{{{\text{st}}}}^{*}\) in Eq. (33) is adopted. For \(\Delta \phi\), Eq. (34) is adopted to avoid direct evaluation of \(K_{{{\text{dyn}}}}^{*}\). After adding and subtracting the static contributions from the sources within the time window, \(\phi _{{\text{w}}}\) can be written as

$$\begin{array}{*{20}c} {\phi _{{\text{w}}} \left( {x,t;t_{{\text{w}}} } \right) = \mathop \sum \limits_{{n = - \infty }}^{\infty } K_{{{\text{st}}}}^{*} \left( n \right)D_{n} \left( t \right){\text{e}}^{{2{{\pi \text{i}}}nx/2L}} } \\ { + \frac{1}{{2{{\pi }}}}\mathop \smallint \limits_{{ - \infty }}^{\infty } {\text{d}}k\mathop \smallint \limits_{{t - t_{{\text{w}}} }}^{t} {\text{d}}t^{\prime}\left[ {\left( {K_{{{\text{dyn}}}} \left( {k,t - t^{\prime}} \right) + K_{{{\text{st}}}} \left( k \right)} \right)\dot{D}\left( {k,t^{\prime}} \right){\text{e}}^{{{\text{i}}kx}} } \right]} \\ { - \mathop \sum \limits_{{n = - \infty }}^{\infty } \mathop \smallint \limits_{{t - t_{{\text{w}}} }}^{t} {\text{d}}t^{\prime}\left[ {K_{{{\text{st}}}}^{*} \left( n \right)\dot{D}_{n} \left( {t^{\prime}} \right){\text{e}}^{{2{{\pi \text{i}}}nx/2L}} } \right]~.} \\ \end{array}$$
(36)

Note that the argument of \(D_{n}\) in the first term has changed from \(t - t_{{\text{w}}}\) in Eq. (33) to \(t\). As mentioned earlier, discretization of the inverse Fourier transform (integral with respect to \(k\)) with a regular mesh causes artificial periodic boundaries every \(\lambda\). If \(\lambda\) is set to a long enough value so that the elastic waves from replication do not reach \(\left[ { - L/2,L/2} \right]\) within the time window \(t_{{\text{w}}}\), then the artificial periodic boundaries do not matter in the second term in Eq. (36). Because \(K_{{{\text{dyn}}}}^{*}\) is not used and the static kernel \(K_{{{\text{st}}}}^{*}\) is very simple in all modes, implementation of this method is rather straightforward for all modes.

The second term in Eq. (36) is the full stress transfer from a source within the time window \(\Delta \phi\) (Eq. (34)). Following Eq. (3), it can be decomposed as \(\Delta \phi = \Delta \phi _{{{\text{dyn}}}} + \Delta \phi _{{st}}\), where \(\Delta \phi _{{{\text{dyn}}}}\) and \(\Delta \phi _{{{\text{st}}}}\) are contributions of \(K_{{{\text{dyn}}}}\) and \(K_{{{\text{st}}}}\). For a compact source, \(\Delta ~\phi\) has a causality cone so that \(\Delta \phi = 0\) outside it. Because the \(\Delta \phi _{{{\text{st}}}}\) reaches infinitely far from the source, \(\Delta \phi _{{{\text{dyn}}}} = - \Delta \phi _{{{\text{st}}}}\) outside the causality cone. If the inverse Fourier transformation is executed with a regular mesh in the wavenumber domain, \(\Delta \phi _{{{\text{dyn}}}} ~\) has the effect of periodic replication outside the causality cone. Therefore, an apparently reasonable formulation

$$\begin{array}{*{20}c} {\phi _{{\text{w}}} \left( {x,t;t_{{\text{w}}} } \right) = \mathop \sum \limits_{{n = - \infty }}^{\infty } K_{{{\text{st}}}}^{*} \left( n \right)D_{n} \left( t \right){\text{e}}^{{2{{\pi \text{i}}}nx/2L}} } \\ { + \frac{1}{{2{{\pi }}}}\mathop \smallint \limits_{{ - \infty }}^{\infty } {\text{d}}k\mathop \smallint \limits_{{t - t_{{\text{w}}} }}^{t} {\text{d}}t^{\prime}\left[ {K_{{{\text{dyn}}}} \left( {k,t - t^{\prime}} \right)\dot{D}\left( {k,t^{\prime}} \right){\text{e}}^{{{\text{i}}kx}} } \right]~} \\ \end{array}$$
(37)

cannot remove the effect of the artificial periodic boundaries.

Numerical implementation

The basic numerical procedure is the same as that used by Noda and Lapusta (2010). A rate- and state-dependent friction law is considered, and the adaptive time-stepper developed by Lapusta et al. (2000) is used. The minimum timestep \(\Delta t_{{\min }}\) is defined during which the S-wave travels a fraction of the grid interval, and the time step is constrained to be an integral multiple of \(\Delta t_{{\min }}\). The time integral in Eq. (36) is calculated using the mid-point rule. Therefore, evaluation of the dynamic kernel is needed only at discrete delay times \(\left\{ {0.5\Delta t_{{\min }} ,~~1.5\Delta t_{{\min }} ,~~2.5\Delta t_{{\min }} ,~~3.5\Delta t_{{\min }} , \ldots ~} \right\}\), which are calculated before the main part of the simulation and stored in memory. A second-order predictor–corrector method is adopted as the time-marching scheme. For details, please see Lapusta et al. (2000) and Noda and Lapusta (2010).

The region of non-zero slip \(\left[ { - L/2,L/2} \right]\) is discretized by using \(N\) regularly spaced grid points at intervals of \(\Delta x = L/N\). The full space \(\left[ { - L,L} \right]\) has \(2N\) grid points. Accordingly, the summation of the infinite series in Eq. (36) is truncated from \(- N + 1\) to \(N\), and the Fourier series \(D_{n}\) and \(\dot{D}_{n}\) are calculated with FFT from arrays of \(\delta ^{\prime}\) and \(V^{\prime}\) of length \(2N\).

To calculate the second term on the right-hand side of Eq. (36) using \(\dot{D}_{n}\), the integration kernel in the spectral domain is prepared in the following way:

  • A slightly longer region than \(\left[ {0,L + ct_{{\text{w}}} } \right]\) is considered and discretized by grid points with the same interval \(\Delta x = L/N\), where \(c = c_{{\text{s}}}\) for an antiplane problem and \(c = c_{{\text{p}}}\) for an inplane problem.

  • A spatio-temporally isolated unit source of \(V^{\prime}\) (in which only one grid point has a unit relative displacement rate for only one time step) is assumed and the corresponding \(\phi\) (\(\phi _{{{\text{aux}}}}\)) is calculated using the spectral method by Lapusta et al. (2000).

  • \(\phi _{{{\text{aux}}}}\) for the relative location from \(0\) to \(\Delta x\left( {N - 1} \right)\) is extracted, and a \(2N\)-long array including the stress change at the relative location from \(- \Delta x\left( {N - 1} \right)\) to \(\Delta x\left( {N - 1} \right)\) is constructed by assuming symmetry.

  • The FFT of the \(2N\)-long array gives the spectral integration kernel in the second term in Eq. (36).

The second and the third terms in Eq. (36) can be calculated at once by subtracting \(K_{{{\text{st}}}}^{*} \left( n \right)\Delta t_{{\min }}\) from the integration kernel for the second term. In the conventional SBIEM, the system size \(\lambda\) must be taken long enough to avoid the effect of periodic replication in \(\Delta \phi\). The extraction of the part of \(\phi _{{{\text{aux}}}}\) that matters in simulation of the fault of length \(L\) enables us to remove the periodic boundary without increasing the system size.

The proposed method has some major advantages compared with the previous SBIEM in earthquake sequence simulation: there is no effect of an artificial periodic boundary condition, and the new method is more efficient. The length of the calculated domain is \(2L\), and the number of grid points in the spatial and wavenumber coordinates is \(2N\). In the simulation by Lapusta et al. (2000), for example, a region with size \(4L\) is used, and thus \(4N\) grid points are considered to reduce the effect of the artificial periodic boundaries. Therefore, with the new method, the memory requirement is reduced by a factor of about 2. Because the time-marching scheme is the same, the smaller number of grid points directly leads to less calculation and shorter CPU time per time step. The integration kernel is calculated before the main part of the simulation, and the additional calculation cost of the preparation of the kernel is not significant.

Verification and evaluation of the effect of periodic boundaries

Problem setting

To verify the proposed method and to evaluate the effect of periodic boundaries, a rather simple test problem is solved with the new method and with the conventional method (Lapusta et al. 2000; Noda and Lapusta 2010), and then the results are compared. In the description of the problem below, scales of stress and speed are given because parameters such as the shear modulus, the S-wave speed, and the plate convergence rate are relatively well constrained for many rocks and major plate boundaries. The length scale can be selected as another independent dimension. Because of the broad size-distribution range of earthquakes, the length scale of interest depends on the specific problem being addressed. Here, \(L = 2~{\text{km}}\) was selected only for presentational purposes. For a longer or shorter fault, the solution does not change if the parameters and variables with the dimension of length are appropriately scaled. For example, for a longer fault \(L = 20~{\text{km}}\), length and time (= length × speed) (e.g., \(x\), \(\delta ^{\prime}\), \(t\)) should be multiplied by 10.

Here, an antiplane problem for an elastic medium with \(\mu = 30~{\text{GPa}}\) and \(c_{\text{s}} = 3000~{\text{m}}/{\text{s}}\) is considered. A sequence of earthquakes on a planar fault with a long-term slip rate of \(V_{{{\text{pl}}}} = 10^{{ - 9}} ~{\text{m}}/{\text{s}}\) is simulated. As the reference, the situation of a uniform slip rate \(V_{{{\text{pl}}}}\) over the entire fault is selected so that \(\tau _{0}\) is independent of \(t\), because a uniform slip on an infinite planar fault causes no stress change. The actual slip rate \(V\) is given by

$$\begin{array}{*{20}c} {V = V_{{{\text{pl}}}} + V^{\prime}~.} \\ \end{array}$$
(38)

The fault is governed by the aging law, a rate- and state-dependent friction law (Dieterich 1979; Ruina 1983):

$$\begin{array}{*{20}c} {\tau = \tau _{0} + A\log \left( {\frac{V}{{V_{{{\text{pl}}}} }}} \right) + B\log \left( {\frac{{V_{{{\text{pl}}}} \theta }}{{d_{{\text{c}}} }}} \right)~,} \\ \end{array}$$
(39)
$$\begin{array}{*{20}c} {\frac{{{\text{d}}\theta }}{{{\text{d}}t}} = 1 - \frac{{V\theta }}{{d_{{\text{c}}} }}~,} \\ \end{array}$$
(40)

where \(A\) and \(B\) represent the amplitudes of the direct effect and evolution effect, respectively, \(\theta\) is the state variable, and \(d_{c}\) is the characteristic slip of the state evolution. The initial value of \(\theta\) is set at \(10d_{c} /V_{{{\text{pl}}}}\). Distributions of \(A\), \(B\), and \(d_{{\text{c}}}\) are plotted in Fig. 2. The mirror symmetry is broken by assuming a gradient of \(B\) to avoid simulation near an unstable symmetric solution; \(B\) is given by multiplying a linear function by an infinitely differentiable smoothed boxcar function (e.g., Noda and Lapusta 2010). The fault is rate weakening \(\left( {A < B} \right)\) in its central part, the length of which is about 1 km. There, \(A = 1.2~{\text{MPa}}\), \(B \approx 1.6~{\text{MPa}}\), and \(d_{c} = 1~{\text{mm}}\), which yields \(A/B \approx 0.75\) and a nucleation size of 191 \(~{\text{m}}\) (Rubin and Ampuero 2005), about one-fifth the length of the rate-weakening region. The length scale of the cohesive zone \(\Lambda _{0}\) (Lapusta and Liu 2009) is about \(16.6~{\text{m}}\), and \({{\Delta }}x = 5~{\text{m}}\) is used in the simulations. The minimum time step and the time window for the inertial effect are set to \(\Delta t_{{\min }} = 0.5\Delta x/c_{{\text{s}}}\) and \(t_{w} = 3L/c_{{\text{s}}}\), respectively.

Fig. 2
figure 2

Distributions of the frictional properties along the fault in the test problem used to verify the proposed method. \(A\) and \(B\) are the amplitudes of the direct effect and evolution effect, respectively, and \(d_{c}\) is the characteristic slip of the state evolution

Full size image

Effects of artificial periodic boundaries

Figure 3 shows the simulation results without periodic boundaries (Fig. 3a) and with periodic boundaries at intervals of \(\lambda = 2L\) (Fig. 3b) and \(\lambda = 10L\) (Fig. 3c). The coseismic period is defined by using \(0.1~{\text{m}}/{\text{s}}\) as the spatially maximum slip rate threshold, and the fault behavior is simulated till the end of the fifth earthquake. In this example, the case with \(\lambda = 10L\) looks quite similar to the case without periodic boundaries (Fig. 3a, c).

Fig. 3
figure 3

Simulation results for a sequence of five earthquakes a without periodic boundaries, b \(\lambda = 2L\), and c \(\lambda = 10L\), where \(L = 2~{\text{km}}\). Blue and red lines indicate distribution of slip every \(10^{7}\) s and \(10^{{ - 2}}\) s during earthquakes, respectively. The case in a was calculated using the proposed method, and those in b and c were calculated by the method of Noda and Lapusta (2010)

Full size image

The case with \(\lambda = 2L\) has larger coseismic slip and thus longer recurrence intervals and a longer simulated time. In addition, the creep fronts penetrate into the rate-weakening patch more significantly. In the third earthquake, dynamic rupture starts not by regular nucleation but by coalescence of creep fronts at the center of the patch, which is known to occur when the brittleness of a rate-weakening patch is low (e.g., Chen and Lapusta 2009). Such irregular nucleation takes place only with \(\lambda = 2L\) as long as studied. Nucleation always occurred in the right part of the rate-weakening patch for \(\lambda \ge 3L\).

As a measure of the effect of the artificial periodic boundaries, the interval between the fourth and the fifth earthquakes is focused on. In the case without periodic boundaries, this interval is \(T_{\infty } = 3.24 \times 10^{8} ~{\text{s}}\). In the cases with periodic boundaries, the interval \(T_{{\lambda /L}}\) is longer than \(T_{\infty }\) and depends on \(\lambda /L\). Figure 4 shows the effect of artificial extension of the recurrence interval as a function of \(\lambda /L\). The normalized extension \(\left( {T_{{\lambda /L}} - T_{\infty } } \right)/T_{\infty }\) decreases with increasing \(\lambda /L\) approximately following \(\left( {\lambda /L} \right)^{{ - 2}}\). This result is consistent with stress perturbation due to a dislocation pair at a distance. The convergence of \(\left( {T_{{\lambda /L}} - T_{\infty } } \right)/T_{\infty }\) towards zero indicates that the proposed method appropriately produces the limit solution of an infinitely large \(\lambda\) obtained by the previous method.

Fig. 4
figure 4

Artificial extension of the earthquake recurrence interval due to periodic boundaries normalized by the case without periodic boundaries. The intervals between the fourth and fifth earthquakes were used. The black points connected by black lines are simulated points, and the dashed line indicates dependency of (λ/L)−2

Full size image

The effect of the artificial periodic boundaries depends on the setting of the problem. Its comprehensive evaluation for different problem settings or different sets of specific parameters is outside the scope of this study. However, the values of the physical properties in the test problem are typical of those used in previous studies (e.g., Lapusta et al. 2000). Therefore, it is possible to obtain a rough sense of the likely effect of using artificial periodic boundaries in those previous studies from Fig. 4. The case without periodic boundaries uses the same amount of computational resources as the case with \(\lambda = 2L\), which yields a recurrence interval longer by 3.6%. The case with \(\lambda = 4L\), which was adopted by Lapusta et al. (2000), overestimates \(T_{\infty }\) by a little less than 1% (0.71%). To decrease the overestimation to less than 1‰, \(\lambda > 10L\) is required, but that would cost significantly more than the proposed method. The longer recurrence interval for shorter \(\lambda\) is likely caused by smaller loading rate to the seismogenic fault patch due to more dense periodic replications of the locked patch. Even with \(\lambda = 2L\), the overestimation is not significant compared with the variability and uncertainty of the recurrence intervals of large earthquakes. Therefore, the use of artificial periodic boundaries would probably have a negligible effect on the conclusions of previous studies.

Discussion and conclusion

In this paper, a new SBIEM method for simulation of dynamic earthquake sequence that is free from the periodic boundary condition is proposed. The new method differs from the conventional method (e.g., Lapusta et al. 2000) with respect to system size and the integration kernel.

With the conventional method, the artificial periodic boundary interval \(\lambda\) is taken to be long (e.g., \(4L\) in Lapusta et al. 2000) to reduce the effect of periodic replication. For a potential non-zero source region of length \(L\), the new method requires a system size no longer than \(2L\). Because this length is shorter than the typical system size in previous studies, it saves computational costs.

In the new method, the full wave-mediated stress transfer in the recent history within the time window \(t^{\prime} \ge t - t_{{\text{w}}}\) is calculated by the conventional method to avoid numerical evaluation of a 2-dimensional integration. The integration kernel is calculated using an auxiliary long domain with realization of a discretized isolated source. Because of the truncation of the temporal integral and the existence of a causality cone in the elastodynamics, the use of a sufficiently long auxiliary domain effectively removes the effect of an artificial periodic boundary. The static effect extends infinitely far from the source and thus requires removal of the artificial periodic boundaries by the technique proposed by Cochard and Rice (1997). This technique causes an integral in the static kernel that can be executed analytically. The new SBIEM is readily applicable to 2-dimensional problems of all rupture modes. As pointed out by Cochard and Rice (1997), the new method can be used even for 3-dimensional problems with a flat fault, because such problems can be formulated by the superposition of 2-dimensional problems (Lapusta and Liu 2009).

To verify the proposed method and to evaluate the effect of artificial periodic boundaries, a simple test problem of an earthquake sequence was solved using both the proposed method and the conventional method with various \(\lambda\). The physical parameters of the problem were set to typical values so that the result would be relevant to many previous models. The artificial periodic boundary increases coseismic slip and the recurrence interval of earthquakes. Even with \(\lambda = 4L\), the recurrence interval is overestimated by only about 1%. This effect is so minor that the conclusions of the previous studies are not affected. In addition, small \(\lambda\) causes long penetration of the creep front before nucleation.

If a problem without a periodic boundary condition is considered, the present method is advantageous over the conventional method not only with respect to the quality of the solution, but also with respect to the computational cost.