2D numerical models of ductile rupture propagation

Earthquakes are commonly associated with brittle failure and frictional sliding in the uppermost 70 km of the Earth. Yet, a significant fraction of seismic events are detected at depths of up to 700 km (deep earthquakes). As these events are difficult to reconcile with our understanding of brittle failure, they are likely facilitated by a ductile weakening mechanism instead. Thermal runaway describes the positive feedback loop of shear heating, temperature-dependent viscosity and deformation. This mechanism has been proposed as a driver of deep earthquakes, and several one-dimensional (1D) studies support its viability. However, two-dimensional (2D) models that show the transient propagation of highly localized shear zones due to thermal runaway are still missing.

We present 2D thermomechanical models which employ a composite visco-elastic rheology, combining elasticity with diffusion creep, dislocation creep and low-temperature plasticity. The code is written in the Julia programming language, operates on Graphic Processing Units (GPU) and utilizes the pseudo-transient relaxation method. Our models capture the nucleation of ductile ruptures on small perturbations, and their transient propagation through a previously intact host rock. Slip velocities inside the ductile ruptures are initially on the order of the far-field deformation, but as the rupture self-localizes, velocities quickly increase by several orders of magnitudes and reach the range of earthquakes (> 1 mm s^-1). The ductile ruptures propagate parallel to the simple-shear background deformation without pre-existing faults or weak layers. If multiple perturbations are present, thermal runaway may nucleate in multiple locations and ruptures can bend to connect to each other.

The magnitude of maximum slip velocity strongly depends on the ratio of stored elastic energy to thermal energy when deformation transitions from low-temperature plasticity to diffusion or dislocation creep. This ratio is derived from one-dimensional models but retains its validity in 2D. If it is small (e.g., low stress, high temperature), shear zones are broad, and deformation is slow. For medium values, slip velocities are in the range of aseismic slow slip events (SSEs). For large energy ratios (e.g., high stress, low temperature), slip velocities reach the seismic window.

Such high-stress conditions are most likely to occur in the cold cores of subducting slabs when they approach the bottom of the mantle transition zone. The resistance of the lower mantle causes slabs to deform, and the large overburden pressure increases viscosity. Both effects increase stress levels. This depth also coincides with the highest occurrence rate of deep-focus earthquakes.