A Divide-and-Conquer Strategy for Fast Elastodynamic Simulation of Earthquakes and Aseismic Slip on Fault Networks

Simulating long-term, fully dynamic sequences of earthquakes and aseismic slip (SEAS) on complex fault networks remains computationally demanding due to the cost of resolving elastodynamic fault interactions. Although high-performance computing environments make such simulations more feasible, the associated expense remains substantial, particularly for multicycle simulations and long-term deformation. Consequently, many studies adopt quasi-dynamic formulations that approximate elastodynamic wave-mediated stress transfer through radiation damping, enabling longer simulations at reduced computational cost. Here, we present an efficient numerical framework for fully elastodynamic SEAS simulations on complex fault networks. The method, presented here in a 2D antiplane setting, follows a divide-and-conquer strategy in which elastodynamic self-effects and fault-to-fault interaction effects are treated separately and evaluated using boundary integral formulations tailored to each interaction type. Self-interactions along planar fault segments are efficiently computed through a nonreplicating spectral boundary integral formulation, eliminating spurious contributions from periodic fault images. Stress transfer across arbitrarily oriented faults is instead resolved using a fully dynamic space–time boundary integral representation accelerated by hierarchical matrices (H matrices). This decomposition defines a general formulation that is readily extensible beyond the present 2D antiplane setting. A central methodological advance is the introduction of a selective H-matrix compression strategy based on the faultwise assembly of independent binary trees, enabling targeted low-rank approximation of long-range fault-to-fault interactions while preserving near-field accuracy and excluding self-interaction effects from the hierarchical structure. Additional computational gains arise from physics-informed truncation of elastodynamic histories, including mode-dependent time-window truncation for spectral self-effects and causality-based truncation for interaction kernels. Benchmark multifault simulations validate the accuracy of the formulation against reference uncompressed space–time solutions. The method reduces the computational and memory complexity of interaction calculations from O(N3) to O(N2logN), yielding up to two orders of magnitude reduction in computation time and approximately one order of magnitude reduction in memory usage for typical problem sizes (∼3·103 degrees of freedom). These gains enable fully dynamic SEAS simulations on complex fault networks using workstation-class hardware.