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

TL;DR

This paper introduces an efficient divide-and-conquer numerical framework for simulating complex earthquake and aseismic slip sequences on fault networks, significantly reducing computational costs and enabling fully dynamic elastodynamic simulations on standard hardware.

Contribution

The authors develop a novel boundary integral method with hierarchical matrix acceleration and physics-informed truncation, improving simulation speed and memory efficiency for complex fault networks.

Findings

Achieves up to three orders of magnitude speedup.

Reduces memory usage by an order of magnitude.

Validates accuracy against reference solutions.

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.