Understanding the physics of thermal runaway and ductile rupture propagation

Ductile deformation is commonly associated with slow and uniform deformation which unfolds over thousands to hundreds of millions of years. Nevertheless, ductile instabilities can result in the localization of deformation into narrow shear zones which operate on much shorter time scales. Mylonites are one example of ductile localization, and events such as slow slip events and deep earthquakes are also associated with fast ductile deformation. The latter are reported up to depths of 700 km and are difficult to reconcile with our understanding of brittle failure which suggests that they are driven by a ductile localization mechanism instead.

One such mechanism is thermal runaway, a feedback loop of shear heating, temperature-dependent viscosity and deformation. Several one-dimensional (1D) studies support the viability of thermal runaway as a driver of deep earthquakes. Here, we present two-dimensional (2D) thermomechanical models of thermal runaway in olivine under simple and pure shear conditions in line with the cold cores of subducting slabs. The models employ a composite visco-elastic rheology including diffusion creep, dislocation creep, and low-temperature plasticity.

The code is written in the Julia programming language and utilizes the package ParallelStencil.jl for GPU parallelization as well as JustPIC.jl for particle-in-cell advection. We employ an accelerated pseudo-transient (APT) solver which makes use of recent developments in automatic tuning of numerical parameters such as pseudo-time steps and damping coefficients. These features enable us to locally employ strong grid refinement (factor 100) without destabilizing the solver which allows us to use model domains spanning tens of kilometers with local grid resolutions of up to one meter.

Our models capture the nucleation and transient propagation of ductile ruptures through a previously intact host rock. The ruptures initiate in zones of reduced grain size and under plate tectonic deformation rates (10 cm yr^-1). During propagation, they self-localize and accelerate to reach slip velocities in the range of earthquakes (> 1 mm s^-1). The magnitude of maximum slip velocity is strongly coupled to the stress in the host rock prior to rupture nucleation, and the ruptures self-consistently run out in the high-temperature/low-stress areas of the model domain. This behavior is consistent with scaling laws derived from 1D models and the occurrence of deep-focus earthquakes in the cold olivine cores of subducting slabs. As we consider the latent heat of melting, our models demonstrate that the local temperature surge due to thermal runaway is fast enough to completely melt a thin layer of olivine during rupture propagation, indicating a link between deep earthquakes and pseudotachylytes.