  Earthquake fault zones are more complex, both geometrically and
rheologically, than an idealised infinitely thin plane embedded in linear
elastic material. To incorporate nonlinear material behaviour, natural
complexities, and multi-physics coupling within and outside of fault zones,
here we present a first-order hyperbolic and thermodynamically compatible
mathematical model for a continuum in a gravitational field which provides a
unified description of nonlinear elasto-plasticity, material damage and of
viscous Newtonian flows with phase transition between solid and liquid phases.
The fault geometry and secondary cracks are described via a scalar function
$\xi \in [0,1]$ that indicates the local level of material damage. The model
also permits the representation of arbitrarily complex geometries via a diffuse
interface approach based on the solid volume fraction function $\alpha \in
[0,1]$. Neither of the two scalar fields $\xi$ and $\alpha$ needs to be
mesh-aligned, allowing thus faults and cracks with complex topology and the use
of adaptive Cartesian meshes (AMR). The model shares common features with
phase-field approaches but substantially extends them. We show a wide range of
numerical applications that are relevant for dynamic earthquake rupture in
fault zones, including the co-seismic generation of secondary off-fault shear
cracks, tensile rock fracture in the Brazilian disc test, as well as a natural
convection problem in molten rock-like material.
