A divide-and-conquer strategy for fast elastodynamic simulation of earthquakes and aseismic slip on fault networks

Abstract

Simulating long-term, fully dynamic sequences of earthquakes and aseismic slip (SEAS) on geometrically complex fault networks remains computationally demanding due to the cost of resolving elastodynamic interactions. Although high-performance computing improves feasibility, simulations remain expensive, particularly for multicycle evolution, motivating the widespread use of quasi-dynamic approximations based on radiation damping. Here we present an efficient numerical framework for fully elastodynamic SEAS simulations on complex fault networks. The method adopts a divide-and-conquer strategy in which elastodynamic self-effects and fault-to-fault interactions are treated separately using boundary integral formulations tailored to each interaction type. Self-interactions along planar faults are computed using a non-replicating spectral boundary integral formulation that eliminates periodic-image artifacts, while interactions between arbitrarily oriented faults are evaluated through a fully dynamic space-time boundary integral representation accelerated by hierarchical matrices (H-matrices). A key advance is a selective H-matrix compression strategy based on fault-wise assembly of independent binary trees, enabling low-rank approximation of long-range interactions while preserving near-field accuracy and excluding self-effects from the hierarchical structure. Additional efficiency arises from physics-informed truncation of elastodynamic histories using mode-dependent time windows and causality-based kernel truncation. Benchmark multi-fault simulations validate accuracy against reference uncompressed solutions. The method reduces interaction complexity from O(N^3) to O(N^2 log N), yielding up to three orders of magnitude speedup and an order-of-magnitude memory reduction for typical problem sizes (~3e10 degrees of freedom), enabling fully dynamic SEAS simulations on workstation hardware.

Figures (15)

• Figure 1: Schematic representation of a fractured domain composed of pre-existing planar faults arbitrarily oriented in space. Fault segments $\Gamma_i$ may be spatially separated or intersect with one another and interact elastodynamically through stress transfer. In the present study, fault slip is governed by Rate-and-State friction and is modeled in an anti-plane (mode III) setting, for which the problem is scalar and only shear tractions are involved.
• Figure 2: Illustration of the hierarchical matrix ($\mathcal{H}$-matrix) construction used to accelerate fault--to--fault interaction calculations. Left: each fault is decomposed into a binary tree of spatial clusters down to leaf nodes. Interactions between well-separated clusters are identified as admissible and represented using low-rank approximations, while near-field interactions are stored in full. Right: schematic representation of the resulting stack of $\mathcal{H}$-matrices associated with the discretized space--time convolution, with one matrix per convolution time lag. Green blocks denote admissible interactions stored in low-rank form, whereas red blocks correspond to near-field interactions stored as full, dense matrices.
• Figure 3: Geometry of the two-fault problem, with all the dimensions expressed in kilometers. Two disconnected planar faults, $\Gamma_1$ and $\Gamma_2$, with arbitrary orientations are embedded in an unbounded elastic medium. Local fault reference systems are defined consistently so that signed slip, slip rate, and shear traction follow the same convention on both faults. Blue profiles indicate the imposed initial shear stress distributions: uniform on fault $\Gamma_2$ and locally peaked on fault $\Gamma_1$ to artificially trigger a dynamic rupture.
• Figure 4: Comparison between the proposed hybrid spectral / space--time formulation with $\mathcal{H}$-matrix acceleration (black solid lines) and the classical uncompressed space--time BIEM (red dashed lines). Results are shown for shear traction, slip rate, and slip at selected time snapshots on (a) fault $\Gamma_1$ and (b) fault $\Gamma_2$. The two solutions are in excellent agreement, indicating no observable loss of accuracy due to hierarchical compression.
• Figure 5: Computational performance comparison for the two--fault problem as a function of the total number of fault elements $N_{\mathrm{elts}} = N_{\mathrm{elts},\Gamma_1} + N_{\mathrm{elts},\Gamma_2}$. (a) Computation time per time step (in seconds). (b) Memory used to store pre-integrated interaction kernel. For sufficiently large $N_{\mathrm{elts}}$, the hybrid spectral / space--time formulation with $\mathcal{H}$-matrix acceleration achieves approximately three orders of magnitude reduction in computation time and one order of magnitude reduction in memory usage compared to the classical space--time BIEM. Pink regions denote the simulation results obtained with cohesive zone length-scales resolved with less than 10 elements.