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.