Why we care: Earthquakes rank among the most destructive natural hazards, capable of claiming thousands of lives and causing devastation to infrastructure and communities within seconds. Understanding when and where large earthquakes are likely to occur is therefore not only a scientific challenge but a pressing societal need. At the heart of this challenge lies the fundamental cyclic nature of earthquakes: tectonic stresses accumulate slowly on faults over decades to centuries, then release abruptly during rupture — only to begin building again. This cyclic character carries a profound implication: if the controlling physical processes and the current stress state of a fault system are sufficiently well constrained, the cycle's repetitive nature bears the potential to inform forecasts of future fault behavior. However, cyclic does not mean periodic. Unlike a clock, fault systems do not rupture on a fixed schedule — recurrence intervals vary considerably even on the same fault, governed by the complex interplay of fault geometry, friction, loading rates, and stress interactions with neighboring faults. This variability is a fundamental property of fault systems that must be understood and accounted for in any serious hazard assessment.
What we do: The non-periodic nature of earthquake recurrence has an important and often underappreciated consequence: even faults with average recurrence intervals of centuries or millennia cannot be considered inactive at any given moment. The timing of the last large rupture is frequently poorly constrained, and even where it is known, the elapsed time since that event tells us relatively little without knowing where the fault currently sits within its stress cycle. A fault may be late in its cycle, or may have experienced accelerated stress loading due to a neighboring event — meaning that hazard assessments based on past earthquake occurrence alone can be misleading. Earthquake cycle simulations address this challenge by modeling the full sequence of interseismic stress accumulation, coseismic rupture, and postseismic relaxation in a physically self-consistent framework. Run over many successive cycles and under varying boundary conditions, these models generate long synthetic earthquake catalogs that far exceed what the observational record alone can provide, allowing us to explore the full range of plausible behaviors a fault system may exhibit — including the rare but inevitable large events that define its hazard.
Our tools: Our primary tool for earthquake cycle simulations is MCQsim (Multi-Cycle earthquake simulator), an open-source simulator developed in our group. MCQsim is written in C and parallelized via MPI, enabling simulations that scale from personal computers to large HPC clusters, and supports fault systems of arbitrary geometric complexity discretized into triangular elements. Efficiency at large scale is achieved through QuadTree and H-matrix compression of the stiffness kernel, which dramatically reduces both memory requirements and computation time — enabling simulations with hundreds of thousands of fault elements that would otherwise be computationally prohibitive. MATLAB-based graphical interfaces handle pre- and post-processing, keeping the simulation workflow accessible without requiring users to work through raw binary files. MCQsim is continuously updated (improving performance and rupture physics representation) and freely available at github.com/OlafZielke-EQ/MCQsim.
Above: Simulated accumulation slip due to M6+ earthquakes along the Gulf of Aqaba fault system (2,000-year time window).
Predicting how often large earthquakes occur on a given fault requires knowing the physical conditions that govern rupture — friction properties, temperature and pressure at depth, fault geometry, and long-term loading rates. In practice, these parameters are rarely well constrained. They vary spatially across a fault surface, may evolve over time, and are largely inaccessible to direct measurement. This uncertainty makes it difficult to move from a general understanding of earthquake recurrence to reliable, fault-specific estimates.
Earthquake cycle simulations offer a systematic way to explore this problem. Rather than seeking a single "best" solution, we use our simulation framework to perform controlled parameter space investigations — varying one or more physical conditions while keeping others fixed — to map out how different assumptions translate into different rupture behaviors and recurrence statistics.
This approach allows us to identify which parameters most strongly control recurrence and which play only a secondary role, providing a clearer picture of where observational effort is best directed. It also quantifies the range of recurrence behavior consistent with a given set of boundary conditions, yielding not a single recurrence estimate but a physically grounded distribution — a more honest and ultimately more useful basis for hazard assessment.
Above: Influence of fault strength distribution (left) on earthquake catalog and large-event recurrence.
Seismic hazard assessment requires knowing not just where earthquakes occur, but how large they are likely to be and how frequently — information that directly informs building codes, insurance frameworks, emergency planning, and infrastructure design. Observational catalogs, while invaluable, are too short to reliably characterize the full magnitude-frequency distribution of large, infrequent events on individual faults — precisely where tail risk is greatest and observational data scarcest.
Earthquake cycle simulations address this gap by generating long synthetic catalogs for naturally occurring fault systems, providing physics-based earthquake rupture forecasts that can serve as input to probabilistic seismic hazard analysis (PSHA). Crucially, running simulations under different physical assumptions allows us to propagate uncertainty in fault parameters directly into the hazard estimate — producing not a single forecast but a range of plausible outcomes that honestly reflects what we do and do not know.
The resulting rupture forecasts are internally consistent — magnitude, recurrence, and rupture extent emerge from the same physical model rather than being specified independently. This consistency makes them a meaningful complement to empirically derived inputs in the PSHA logic tree, particularly for faults where observational constraints are sparse.
Above: Toy model earthquake catalog and resulting hazard map for the East Anatolian fault system, showcasing the integration of simulated catalogs into PSHA engines
Energy resource extraction and subsurface fluid injection alter the stress and pore-pressure conditions at depths where earthquakes originate. Even in regions with little natural seismicity, these pressure changes can reactivate existing faults or fractures, triggering earthquakes that would not otherwise have occurred. Understanding and anticipating this induced seismicity is a critical component of responsible resource extraction and geothermal energy development.
Our simulation framework allows us to investigate how sudden or sustained pressure changes affect earthquake probability — both on large, well-characterized fault structures and within discrete fracture networks (DFNs), which more realistically represent the complex, multi-scale fracture systems typically found at reservoir depths in the absence of major tectonic structures. A key question is how the timing and magnitude of fluid injection interact with where a fault or fracture network currently sits within its stress cycle — a small pressure perturbation near the end of a cycle may trigger a large event that would have occurred anyway, while the same perturbation early in the cycle may produce only minor seismicity.
This work helps to define operational thresholds for safe injection practices and to assess seismic risk in both tectonically active and critically stressed but otherwise dormant regions — contributing to the broader goal of developing subsurface energy resources without unacceptable seismic consequences.
Above: Rupture slip for a 2,000-year time window across a discrete fracture network consisting of 1000 individual fractures.
Above: Simulation of induced earthquake activity along two two faults, caused by fluid injection and associated pore pressure changes (after C. Salah, MSc student at KAUST).