How to break the lithosphere: a compressible pressure-velocity formulation for elasto-visco-plastic rheologies that includes shear and tensile failure with dilation

One of the great challenges involved in modelling the lithosphere is its plastic behaviour, especially when dealing with compressible materials. Shear fractures are designated as mode 2 and 3 and can be described using a Linear Mohr Coulomb envelope  or a simplification of it like Drucker-Prager. Meanwhile, mode 1 fractures are created when the normal stresses become tensile  and require another yield function, such as the Griffith criterion or a tension cap function.

While the governing equations are well known and widely employed in engineering codes, they are usually expressed with a displacement formulation. Most geodynamic codes, on the other hand, use pressure and velocity as their primary variables. A numerically robust method that takes all plasticity modes into account in a staggered finite difference discretization remains an open task. Here we present a composite yield function implemented with pressure-velocity formulation, capable of producing produce shear and tensile failure.

We have implemented this in a new code that employed PETSc through the recently updated PETSc.jl Julia interface, while utilizing the automatic differentiation tools in julia. We found this workflow to significantly reduce the development time of complex nonlinear coupled codes.  

We will describe the implementation, propose regularization schemes and discuss benchmark cases and simple applications. We demonstrate Newton convergence for most cases and will discuss different methods to combine multiple plastic flow laws.