Numerical Simulation on Faulting: Microscopic evolution, macroscopic interaction and rupture process of earthquakes

We review the recent researches of numerical simulations on faulting, which are interpreted in this paper as the evolution of the state of the fault plane and the evolution of fault structure. The theme includes the fault constitutive (friction) law, the properties of the gauge particles, the initial phase of the rupture, the dynamic rupture process, the interaction of the fault segments, the fault zone dynamics, and so on. Many numerical methods have been developed: boundary integral equation methods (BIEM), finite difference methods (FDM), finite or spectral element methods (FEM, SEM) as well as distinct element methods (DEM), discrete element methods (again DEM) or lattice solid models (LSM). The fault dynamics should be solved as a complex non-linear system, which shows multiple hierarchical structures on its property and behavior. The researches have progressively advanced since the 1990’s both numerically and physically thanks to high performance computing environments. The interaction at small scales is modeled to provide a large scale property of the fault. The dynamic rupture has been actively studied especially for the effect on the fault geometry evolution or due to the existed fault structure. The (quasi-)static and the initial processes of the fault movement have been also explored in a seismic cycle. The effect of fluid or heat has been also taken into account in the mechanics. All these efforts help us to understand the phenomena and the unified understanding (simulation) over different spacio-temporal scales is more and more expected.

💡 Research Summary

**
This paper provides a comprehensive review of recent advances in numerical simulation of faulting, emphasizing the simultaneous evolution of the fault‑plane state (microscopic) and the fault structure (macroscopic). The authors begin by outlining the fundamental challenge: earthquakes are governed by a highly nonlinear, hierarchical system in which small‑scale processes such as grain‑scale friction, contact mechanics, and fluid‑induced weakening interact with large‑scale features like fault geometry, segmentation, and stress heterogeneity. To capture this complexity, a wide spectrum of numerical techniques has been developed since the 1990s.

Boundary Integral Equation Methods (BIEM) were the first to enable efficient modeling of rupture propagation in an infinite elastic medium, but their ability to incorporate fully nonlinear constitutive laws is limited. Finite Difference Methods (FDM) introduced explicit time stepping for nonlinear friction and allowed adaptive grid refinement, improving resolution of near‑fault gradients. Finite Element (FEM) and Spectral Element Methods (SEM) later provided the flexibility to embed complex material behavior—rate‑and‑state friction, temperature‑dependent weakening, poroelastic coupling—within realistic three‑dimensional fault geometries. FEM excels at handling heterogeneous elasticity, plasticity, and heat conduction, while SEM offers high‑order accuracy for wave‑field calculations, making it ideal for studying high‑frequency radiation during dynamic rupture.

Particle‑based approaches, notably the Discrete Element Method (DEM) and Lattice Solid Models (LSM), bring the microscopic perspective to the fore. By representing the fault zone as an assembly of interacting particles or lattice nodes, these methods resolve nucleation of microcracks, grain crushing, and the evolution of contact friction directly. The resulting micro‑mechanics can be statistically up‑scaled or coupled through multiscale frameworks to inform the macroscopic fault strength and rupture speed. This hierarchical linking is a central theme of the review: small‑scale physics is not merely a boundary condition but a driver of large‑scale fault behavior.

High‑performance computing (HPC) has been a catalyst for progress. The authors describe how GPU acceleration, MPI‑based domain decomposition, and adaptive mesh refinement (AMR) have pushed simulations from a few hundred thousand degrees of freedom to tens of millions, enabling fully dynamic rupture calculations that resolve both the initial nucleation phase and the subsequent propagation across complex fault networks. These capabilities have revealed several key insights. First, the geometry of fault segments—bends, step‑overs, and intersecting strands—can either accelerate or arrest rupture, depending on the orientation of the stress field and the local frictional parameters. Second, the inclusion of fluid pressure and thermal effects produces rapid frictional weakening (so‑called “hydro‑thermal runaway”), which can trigger cascading failures and explain aftershock sequences that are otherwise difficult to reproduce with purely elastic models. Third, the interaction between the quasi‑static loading phase of the seismic cycle and the dynamic rupture phase is now being simulated in a unified framework, allowing researchers to track stress accumulation, sudden slip, and post‑seismic relaxation within a single model run.

The review concludes by identifying open challenges and future directions. A major priority is the development of data‑driven sub‑grid models—leveraging machine learning to translate detailed DEM/LSM results into effective constitutive laws for FEM/SEM simulations. Another is the systematic validation of numerical outputs against high‑resolution geodetic, seismic, and InSAR observations, which will tighten the link between theory and real‑world earthquake behavior. Finally, the authors advocate for fully coupled thermo‑hydro‑mechanical–chemical models that can simulate long‑term fault healing, mineral precipitation, and the evolution of fault zone material properties over multiple seismic cycles. By integrating these advances, the community moves closer to a unified, multiscale simulation platform capable of reproducing the full spatio‑temporal spectrum of earthquake phenomena.