Framework

To solve the earthquake dynamic simulation problem using boundary integral equations, we assume the fault plane is embedded in an elastic half-space or full-space with homogeneous and constant elastic moduli. A constant tectonic loading rate is imposed across the entire fault interface. The elastic stress transfer due to slip on the fault is described in Equations (1) (shear stress) and (2) (normal stress), with the radiation damping assumption [1].

(1)\[\tau_{i} = \tau_{0} - \sum_{j=1}^{n} k_{ij}^s \left( u_{j} - V_{pl} t \right) - \frac{\mu}{2c_s} \frac{\partial u_{i}}{\partial t}\]

(2)\[\bar{\sigma}_{i} = \bar{\sigma}_{0} + \sum_{j=1}^{n} k_{ij}^N \left(u_{j} - V_{pl} t\right)\]

where \(V_{pl}\) is the imposed tectonic slip rate, \(\mu\) is the shear modulus, \(c_s\) is the shear wave speed, and \(u_j\) is the slip at the \(j\)-th element. The kernels \(k_{ij}^s\) and \(k_{ij}^N\) represent the shear and normal stiffness matrices, respectively. The last term in Equation (1) captures radiation damping and approximates inertial effects, adopted to avoid unbounded slip velocity that would otherwise develop in quasi-static models.

To compute \(k_{ij}^s\) and \(k_{ij}^N\), analytical formulas for static stress induced by triangular dislocations in a homogeneous elastic full-space and half-space are employed [2]. Since the goal is to simulate three-dimensional complex non-planar fault geometries, optimizations specific to planar faults (e.g. Fourier-domain construction) are not applicable. Instead, CPU-based multiprocessing or MPI are used to accelerate the computation of Green’s functions.

—

Time-Differentiated Formulation

By differentiating Equations (1) and (2), and considering external loading, we obtain:

(3)\[\frac{d\tau_{i}}{dt} = -\sum_{j=1}^{N} k_{ij}^s \left(V_{j} - V_{pl}\right) + \dot{\tau}_{i} - \frac{\mu}{2c_s} \frac{dV_{i}}{dt}\]

(4)\[\frac{d\sigma_{i}}{dt} = \sum_{j=1}^{N} k_{ij}^N V_{j} + \dot{\sigma}_{i}\]

where \(\dot{\tau}_i\) and \(\dot{\sigma}_i\) are tectonic loading rates.

—

Friction Law

To close the system, we apply the laboratory-derived rate-and-state friction (RSF) law [3, 4]. The friction coefficient under the regularized aging law is:

(5)\[f(V, \theta) = a \sinh^{-1} \left[ \frac{V}{2V_0} \exp\left(\frac{f_0 + b\ln\left(\frac{V_0 \theta}{d_c}\right)}{a}\right) \right]\]

(6)\[\frac{d\theta}{dt} = 1 - \frac{V\theta}{d_c}\]

with \(d_c\) the characteristic slip distance, \(V_0\) the reference slip rate, \(f_0\) the reference friction coefficient, and \(a\), \(b\) the RSF parameters.

We define a transformed state variable:

(7)\[\psi = f_0 + b \ln \left( \frac{V_0 \theta}{d_c} \right)\]

giving the transformed law:

(8)\[\frac{\tau_i}{\sigma_i} = a \arcsin\left( \frac{V_i}{2V_0} \exp\left(\frac{\psi_i}{a}\right) \right)\]

(9)\[\frac{d\psi_i}{dt} = \frac{b}{d_c} \left[ V_0 \exp\left( \frac{f_0 - \psi_i}{b} \right) - V_i \right]\]

—

System of Equations

By substitution, Equations (3), (4), and (9) yield a coupled system of ODEs of dimension \(4N\):

(10)\[\frac{dy}{dt} = f(y)\]

with:

(11)\[y = (\psi_1, \dots, \psi_N, \tau_{1,1}, \dots, \tau_{1,N}, \tau_{2,1}, \dots, \tau_{2,N}, \sigma_1, \dots, \sigma_N)\]

This system is solved using the Dormand–Prince 5th-order Runge–Kutta method with adaptive time stepping [5].

—

Hierarchical Matrix Compression and MPI

Following Börm [6], PyQuake3D implements H-matrix compression in Hmatrix.py with MPI-based parallel acceleration. The H-matrix framework decomposes stiffness matrices into low-rank far-field blocks and dense near-field blocks, organized into cluster and block trees.

Cluster tree construction is based on geometric splitting, while block trees pair clusters to determine admissibility for low-rank approximation. This design enables efficient compression and distributed memory scalability.

—

References

[1]

James R Rice. Spatio-temporal complexity of slip on a fault. Journal of Geophysical Research: Solid Earth, 98(B6):9885–9907, 1993.

[2]

Mehdi Nikkhoo and Thomas R Walter. Triangular dislocation: an analytical, artefact-free solution. Geophysical Journal International, 201(2):1119–1141, 2015.

[3]

James H Dieterich. Modeling of rock friction: 2. simulation of preseismic slip. Journal of Geophysical Research: Solid Earth, 84(B5):2169–2175, 1979.

[4]

Andy Ruina. Slip instability and state variable friction laws. Journal of Geophysical Research: Solid Earth, 88(B12):10359–10370, 1983.

[5]

William H. Press, Saul A. Teukolsky, William T. Vetterling, and Brian P. Flannery. Numerical Recipes: The Art of Scientific Computing. Cambridge University Press, Cambridge, 3rd edition, 2007.

[6]

Steffen Börm, Lars Grasedyck, and Wolfgang Hackbusch. Introduction to hierarchical matrices with applications. Engineering analysis with boundary elements, 27(5):405–422, 2003.