Bridging the Gap Between Millions of Years and Milliseconds: Modeling Earthquake Sequences, Slow Slip, and Splay Fault Rupture in Subduction Zones

Earthquakes and tsunamis occur on a timescale of seconds and are experienced by humans as sudden devastating disasters. However, the tectonic systems that determine where they occur are shaped over millions of years. Deformation in subduction zones is characterized by visco-elasto-plastic interactions between the accretionary prism featuring splay faults, subducting and overriding plate, asthenosphere, and free surface. To understand the present-day seismicity, earthquake cycle, and splay faulting in particular, these deformation processes need to be considered across all time scales. However, numerical models have not been able to resolve the dynamics across both tectonic and earthquake time scales.

We present a novel numerical modeling technique that simulates fully dynamic earthquake sequences and slow slip events in a subduction zone described by a visco-elasto-plastic rheology. Faults form and evolve spontaneously according to heterogeneous, temperature-dependent material parameters and the local stress field during both the initial 4 million years of subduction and the subsequent seismic phase. We employ an invariant formulation of rate- and state-dependent friction and adaptive time stepping to fully resolve all phases of the seismic cycle.

We generate events covering the slip spectrum from aseismic creep to earthquakes with slip rates in the order of m/s and tens of meters of slip. We find that events are largely characteristic despite the potential for deviating rupture paths in the subduction channel. We find that splay faults need to be sufficiently weak to be activated during a megathrust earthquake, since they cannot accumulate stress over time because velocity-strengthening afterslip relaxes their stresses. Dynamic triggering of a splay fault can lead to an early arrest of the megathrust rupture. Such short-term effects alter the long-term deformation compared to a purely geodynamic model by increasing the importance of one splay fault over others. We also observe that trapped seismic waves significantly change the slip distribution in a similar manner as has been found using a dynamic rupture model.

We conclude that our model successfully combines aspects of established geodynamic models and dynamic rupture models, providing a missing link between the long-term and the short-term. When applying this modeling approach to a complex continental setting, the interaction of multiple faults results in further complexities such as clustering. This highlights the potential and versatility of the method for a wide range of tectonic settings.