# Volumetric (dilatant) plasticity in geodynamic models and implications on thermal dissipation and strain localization

Ekeabino Momoh, Harsha S Bhat, Stephen Tait, Muriel Gerbault

► To cite this version:

Ekeabino Momoh, Harsha S Bhat, Stephen Tait, Muriel Gerbault. Volumetric (dilatant) plasticity in geodynamic models and implications on thermal dissipation and strain localization. Geophysical Journal International, 2024, 240 (3), pp.1551-1578. 10.1093/gji/ggae463. hal-05018360

HAL Id: hal-05018360

https://hal.science/hal-05018360v1

Submitted on 3 Apr 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Volumetric (dilatant) plasticity in geodynamic models and implications on thermal dissipation and strain localization

Ekeabino Momoh<sup>1,2,3</sup> Harsha S. Bhat,<sup>2</sup> Stephen Tait<sup>1,3</sup> and Muriel Gerbault<sup>1</sup>

Accepted 2024 December 23. Received 2024 December 19; in original form 2024 April 2

## SUMMARY

Here, we present a new thermomechanical geodynamic, numerical implementation that incorporates Maxwell viscoelastic rheology accounting for temperature-dependent power-law dislocation creep and pressure-sensitive, non-associated Drucker–Prager brittle failure, as well as for volumetric stresses and strains during viscoplastic flow, a departure from the traditional incompressible assumptions. In solving for energy conservation, we incorporate the heat source term resulting from irreversible mechanical deformations, which embodies viscoelastic and viscoplastic work, and by considering the total stress tensor and total inelastic strain rate tensors, including dilatant plasticity effects for lithospheric-scale applications, instead of only the shear terms as is usually assumed for incompressible materials. This form of the work term thus allows to consider, volumetric deformation and to couple the energy equation to the constitutive description, and hence the stress balance, via the evolving temperature field. Code design enables us to switch individual features of this general rheology ‘on or off’ and thus to benchmark this implementation with published numerical experiments of crustal-scale shortening experiments. We investigate whether ‘brittle-plastic’ compressibility can promote or inhibit localization of deformation and thermal evolution during compression for crustal, and upper mantle rheology. For both crustal-scale and lithospheric-scale experiments, we establish that the feedback from volumetric dissipation, while contributing to temperature increase along with shear dissipation, can potentially slow down heat production per unit time, depending on the choice of boundary conditions. Our new implementation can be used to address buckling problems and collision tectonics.

**Key words:** numerical modelling; heat generation and transport; Mechanics, theory, and modelling.

## 1 INTRODUCTION

Strain localization, that is, the concentration of internal deformation along narrow paths of finite thickness in a material during loading, is a key ingredient of plate tectonics (Jacquey et al. 2021). It is typically interpreted as a precursor to material failure (Besson et al. 2010) and characterizes fault zones, ranging in length from a few metres to hundreds of kilometres and up to plate boundary scale, whose activity may persist for millions of years. In the context of plate convergence in geodynamics, both initiation of subduction interfaces and/or lithospheric-scale structures in continental collision zones are example cases where strain localization and the evolution of deformation are central problems to study (Regenauer-Lieb & Yuen 1998; Toth & Gurnis 1998; Stern 2004; Stern & Gerya 2018; Lallemand & Arcay 2021). An underlying issue is that experimentally determined strengths of rocks (primarily mantle peridotites)

are too great for the resistance of a homogeneous lithosphere to be overcome by the stresses generated by buoyant mantle convection, making it difficult to break (McKenzie 1977; Cloetingh et al. 1989; Mueller, & Phillips 1991). Strain localization in the lithosphere is associated to an irreversible process and to the formation of shear bands, which may be compared to real fault zones. The onset of localized deformation in numerical or analogue models is usually controlled either by including some heterogeneity in initial conditions, for example in material properties, and/or some rheologic characteristic that will lead to positive feedback. Thus, the existence of prescribed weak zones such as lithospheric-scale faults inherited from previous tectonic episodes is expected to be relevant (Cloetingh et al. 1982; Gurnis 1992; Toth & Gurnis 1998; Nikolaeva et al. 2010; Baes et al. 2011). Rheologic features prone to causing a positive feedback once localization nucleates such as grain-size reduction (Thielmann et al. 2015), structural softening

<footnote n="1"><sup>1</sup>Geosciences Environnement Toulouse – GET, Université Toulouse III – Paul Sabatier, CNRS UMR 5563, IRD, CNES, F-31400 Toulouse, France.
E-mail: ekeabino.momoh@get.omp.eu</footnote>

<footnote n="2"><sup>2</sup>Laboratoire de Geologie, École Normale Supérieure, CNRS UMR 8538, PSL Research University, Paris, F-75005 Paris, France</footnote>

<footnote n="3"><sup>3</sup>Laboratoire Dynamiques des Fluides Géologiques, Institut de Physique du Globe de Paris, UMR 7154-CNRS, Université de Paris, F-75005 Paris, France</footnote>
(Le Pourhiet 2013; Duretz et al. 2016), thermally activated softening due to shear heating during irreversible deformation (Crameri & Kaus 2010; Thielmann & Kaus 2012; Duretz et al. 2015; Willis et al. 2019), density changes or age offsets between plates (Leng & Gurnis 2011, 2015; Zhang et al. 2021; Zhou & Wada 2021) material anisotropy (Pardoen et al. 2015) or damage rheology (Karrech et al. 2011) have been looked at. Thielmann & Kaus (2012), for example, investigated thermal softening due to shear heating during irreversible deformation and explicitly did not include a pre-existing mechanically weak zone to focus on this thermally activated feedback. However, their system was given an initial thermal heterogeneity with an age offset to facilitate nucleation of localization. Some of these various effects might indeed occur concomitantly in natural situations. We note here that strain localization can also arise naturally from non-associative brittle plasticity (Vermeer & de Borst 1984; Gerbault et al. 1998), as a result of kinematic stress rotation effects.

Geodynamic modelling has sometimes been approached from a fluid-mechanical perspective and assumptions of incompressibility, either elastic or plastic (Gerya et al. 2004; Babeyko & Sobolev 2008; Thielmann & Kaus 2012; Schmalholz et al. 2014; Ruh et al. 2015; Vogt et al. 2017), which means that elastic compressibility and/or dilation during plastic flow are unaccounted for. In studies of lithospheric-scale deformation, plastic compressibility is mostly not taken into account (Kaus 2010; Duretz et al. 2014; Bessat et al. 2020), and elasticity is often neglected in studies over long timescales (Garel et al. 2014; Holt et al. 2015; Patocka et al. 2019; Pajang et al. 2021). However, a number of studies do account for elastic and plastic compressibilities into the constitutive laws for large-scale lithospheric models (Hassani et al. 1997; Gerbault 2000; Gurnis et al. 2004; Burov 2011; Burov et al. 2014; Duretz et al. 2020, 2021; Jacquey & Cacace 2020). From a solid mechanics perspective, rocks are known to accommodate elastic behaviour and undergo volumetric strain during plastic deformation and they are therefore not strictly plastically incompressible (Bridgman 1964; Cook 1970; Alejano & Alonso 2005; Zhao & Cai 2010). At least in the colder parts of the lithosphere, the role of elasticity in stress storage and hence its potential influence on the overall deformational behaviour might actually have been underestimated.

Whether or not elasticity or plasticity are included in constitutive laws for geodynamic modelling, the influence of volumetric (dilatant) plastic deformation as an energy source in the energy equation has not been considered yet by the geodynamic community. It remains an open question whether or not the plastic compressibility of rocks under lithospheric conditions may contribute to localization both from a solid-mechanical and a thermal perspective.

In the current contribution, we present and benchmark a new solid-mechanical numerical implementation that explicitly includes volumetric (dilatant plastic) contributions in the momentum and energy balances. We include specifically a source term in the energy equation that is due to irreversible deformational work. We explore the effect of thermal softening as an agent of localization, but in addition to considering the impact of heating due to irreversible shear deformation, we aim to develop an understanding of the role of volumetric (dilatant plastic) deformation. In practical terms, we compare the effect of including volumetric plastic deformations in the thermal feedback. We compute irreversible deformation and heating with the full stress and stain-rate tensors, or with only shear terms contributing thermodynamically as has been mostly done before (Regenauer-Lieb et al. 2001; Kaus & Podladchikov 2006; Nikolaeva et al. 2010; Thielmann et al. 2015), that is, we can switch on or off the different effects, so that we go beyond

previous work. For convenience and to avoid confusion with general volumetric strains associated with thermoelastic effects, and finite deformations, we refer to heating arising from contributions due to irreversible dilatant plastic strains and volumetric stresses as dilatant heating.

We begin by stating the mechanical and energy conservation laws and then give a step-by-step description of the constitutive laws and solution methods we have developed. We then proceed to benchmark our new implementation against published models by presenting comparative tests assuming viscoelastic, and viscoelastic-viscoplastic rheologies, and compare the relative contributions of shear heating and shear combined with dilatant heating. Furthermore, we apply our approach to the geodynamic context of lithospheric-scale strain localization. Finally, we discuss the implications of this study for future work.

# 2 THEORETICAL FRAMEWORK

## 2.1 Conservation laws

Our starting point is the balance of linear momentum for a compressible viscoelastic-viscoplastic solid medium given by Lemaitre & Chaboche (1994) and Zienkiewicz & Taylor (2005):

$$ \frac{\partial \sigma_{ij}}{\partial x_j} + f_i = \rho \frac{\partial^2 u_i}{\partial t^2}, \qquad (1) $$

where $\sigma_{ij}$ represents components of the Cauchy stress tensor, $f_i$ represent internal body forces: in this case the lithostatic stress state due to the weight of the material given by $\rho g$, with $\rho$ as density, $g$ as gravitational acceleration, $x_j$ represents spatial variables in the Cartesian coordinates; $u_i$ are the components of displacement. A quasi-static approach is used, which means that we account for small accelerations in the system. The mass balance in Lagrangian framework is given by:

$$ \rho_R = \det \left( \delta_{ij} + \frac{\partial u_i}{\partial x_j} \right) \rho, \qquad (2) $$

where $\rho_R$ is the mass density at a point in a reference (undeformed) configuration and $\rho$ is the density at that point in the current (deformed) configuration. The quantity in parenthesis maps a quantity from a reference (undeformed) state to a current configuration (Zienkiewicz & Taylor 2005). Under infinitesimal strain theory (where displacements and strains are infinitesimal), mass balance translates to fixed density. Mechanical dissipation arises from irreversible thermomechanical work, that is, due to viscoelastic and viscoplastic deformation. We can account for heat generation and diffusion for the case of constant thermal conductivity, $\lambda$, in the conservation of energy which reads:

$$ \frac{\partial T}{\partial t} = \alpha_{th} \nabla^2 T + \beta \frac{\sigma_{ij}(\dot{\epsilon}_{ij}^{v} + \dot{\epsilon}_{ij}^{vp})}{\rho C_p}, \qquad (3) $$

where $C_p$ (J kg⁻¹ K⁻¹) is the specific heat capacity and $\alpha_{th}$ (m² s⁻¹) is the thermal diffusivity with $\alpha_{th} = \lambda / \rho C_p$. $\dot{\epsilon}_{ij}^v$ and $\dot{\epsilon}_{ij}^{vp}$ strain rate components from (viscous) creep and viscoplastic deformations, respectively; $\sigma_{ij}(\dot{\epsilon}_{ij}^v + \dot{\epsilon}_{ij}^{vp})$ is the contribution to heating from deformational work (Rittel 1999; Ravichandran et al. 2002) with $\beta$, the Taylor-Quinney coefficient which quantifies the proportion of deformational work which is dissipated into heat. In this paper, we choose a value of 1, that is, all mechanical work is dissipated as heat. Which deformation strain rate contributes to the thermal dissipation depends on which of them is dominant at a given point in
space or time and is not imposed a priori. In addition, thermoelastic effects including thermal expansion and adiabatic heating, and density variation with temperature and pressure are not considered here.

## 2.2 Deformational framework

Our formulation relies on a small (infinitesimal) strain theory which allows an additive decomposition of any contributing strain rate components (Mase 1970; Lemaitre & Chaboche 1994; Zienkiewicz & Taylor 2005; Bower 2009; de Souza Neto et al. 2011; de Borst et al. 2012). Our constitutive law incorporates elastic strain, ductile creep and plastic flow laws with inelastic work being dissipated into heat.

We assume a Maxwell viscoelastic rheology and that the strain rate is a combination of an elastic contribution and a viscous contribution (Turcotte & Schubert 2002; de Borst et al. 2012; Gerya 2019), and Drucker-Prager viscoplasticity. Each time step starts with an elastic trial state that is then corrected via viscous creep activated at non-zero stresses, and followed by plasticity through a yield stress criterion. Such rheology has been widely used in earlier geodynamic models (Hassani et al. 1997; Gerbault 2000; Gerbault et al. 2002; Burov & Cloetingh 2010). The standard approach for implementing viscoelasticity in the geodynamic community is a combination of shear modulus G, viscosity η and time step Δt. Recent considerations involve a correction of the shear modulus to obtain a viscoelastic shear modulus (G<sup>ve</sup>) that depends on the dynamic shear viscosity (η), the elastic shear modulus (G) and the algorithmic time step (Δt), (Duretz et al. 2018, 2021).

Based on small-strain assumption, we compute the stress evolution, creep and plastic flow relying on an additive decomposition of the total strain rate tensor into elastic, viscous and viscoplastic strain rates, represented, respectively, by superscripts ‘e’, ‘v’ and ‘vp’ below:

$$ \dot{\epsilon}_{ij} = \dot{\epsilon}_{ij}^{e} + \dot{\epsilon}_{ij}^{v} + \dot{\epsilon}_{ij}^{vp}, \quad (4) $$

$$ \dot{\epsilon}_{ij} = \dot{\epsilon}_{ij}^{e} + \dot{\gamma}^{v} \frac{\partial \Phi_{F}^{v}}{\partial s_{ij}} + \dot{\gamma}^{vp} \frac{\partial \Phi_{F}^{vp}}{\partial \sigma_{ij}}. \quad (5) $$

$\Phi_F^v$ and $\Phi_F^{vp}$ are respectively viscous (creep) and viscoplastic flow potentials, the derivatives of the viscous and viscoplastic flow potentials with respect to deviatoric stress and total stress components are indicative of viscous and viscoplastic flow directions, respectively; $\dot{\gamma}^v$ and $\dot{\gamma}^{vp}$ are viscous and viscoplastic multipliers in rate form; $s_{ij}$ represent the components of the deviatoric stress tensor with $s_{ij} = \sigma_{ij} - \sigma_{kk}/3$.

### 2.2.1 Elastic rheology

The elastic part of the material response is treated as an isotropic solid characterized by its Young’s modulus, E, and Poisson’s ratio, v. The classical elastic stress strain relation is given by Hooke’s Law (Zienkiewicz & Taylor 2005; Bower 2009):

$$ \sigma_{ij}^{e} = C_{ijkl}^{e} \varepsilon_{kl}^{e}, \quad (6) $$

where $C_{ijkl}^e$ is a fourth-order tensor representing the components of the elastic stiffness tensor of a material, incorporating the elastic constants listed above. $\sigma_{ij}^e$ and $\varepsilon_{kl}^e$ are respectively, components of elastic stresses and strain tensors.

### 2.2.2 Creep (viscous) rheology

It is commonly assumed that, at high temperatures and low stresses, rocks deform viscously without stresses necessarily overcoming a yield criterion, that is, yield criterion is effectively zero (de Souza Neto et al. 2011). The limiting behaviour for the elastic rheology in our constitutive law is therefore high temperature creep. We assume that the dominant creep mechanism is dislocation creep and formulate the viscous deformation as follows:

$$ \dot{\epsilon}_{ij}^{v} = \dot{\gamma}^{v}(\sigma_{ij}^{e}, T) \frac{\partial \Phi_{F}^{v}}{\partial s_{ij}^{e}}, \quad (7) $$

with $\dot{\gamma}^v$, a non-negative quantity specifying a magnitude of viscous flow in rate form and $\partial\Phi_F^v/\partial s_{ij}^e$ specifying the direction of viscous flow taking account of only deviatoric stresses. $\Phi_F^v = \sqrt{J_{II}^e}$ with $J_{II}^e$ representing the invariant of the deviatoric stress tensor, is the viscous flow potential whose direction is given by:

$$ \frac{\partial \Phi_F^v}{\partial s_{ij}^e} = \frac{s_{ij}^e}{2\sqrt{J_{II}^e(s_{ij}^e)}}. \quad (8) $$

The magnitude of viscous flow $\dot{\gamma}^v$ assumes various functional forms depending on the specific problem and material (de Souza Neto et al. 2011). We utilize the power-law functional form of (Boyle & Spence 1983; Rosakis et al. 2000; Poulet & Veveakis 2016):

$$ \dot{\gamma}^{v}(\sigma_{ij}^{e}, T) = (\Phi_{F}^{v})^{m} f(T). \quad (9) $$

$m$ represents the power-law exponent that describes the sensitivity to stress during viscous flow. The temperature dependence is given in the form of (Skrzypek 1993; Kohlstedt et al. 1995; Ranalli 1995):

$$ f(T) = A e^{-\frac{E_a}{RT}}, \quad (10) $$

Here $E_a$ is the activation energy, $R$ is the molecular gas constant, $T$ is the absolute temperature (in Kelvin, K) and $A$ is the pre-exponential factor.

Eq. (10) indicates that viscous creep becomes significant at high temperatures; otherwise, the viscous correction at low temperatures is so small that the effect of elasticity is preserved. Because the viscous correction term is insignificant at low temperatures and low stresses, the feedback to viscous strain rate and deviatoric stresses is small, and the contribution to thermal dissipation evolves similarly (thermoelastic effects are unaccounted for). Therefore, the partitioning between elastic behaviour and ductile creep depends essentially on the temperature. This strategy can be seen as an elastic predictor, followed by a viscous correction step (Jacquey & Cacace 2020).

### 2.2.3 Viscoplastic rheology

The elastic and viscous behaviour or viscoelastic rheology is bounded by a pressure-sensitive frictional plastic yield criterion which has been used to model brittle deformation on a lithospheric scale (Moresi et al. 2007; Babeyko & Sobolev 2008; Burov & Cloetingh 2010; Kaus 2010; Baes et al. 2011; Jacquey & Cacace 2020). Here, we use a Drucker-Prager criterion given by (Drucker & Prager 1952; Alejano & Bobet 2012):

$$ \Phi_{Y}^{DP} = \sqrt{J_{II}^{v}(s_{ij}^{v})} + \alpha_{1} P^{e} - \alpha_{2} c > 0, \quad (11) $$

where $J_{II}^v = s_{ij}^v s_{ij}^v / 2$ represents the second invariant of the deviatoric stress tensor after creep deformation, $P^e = I_1 / 3 = \sigma_{kk} / 3$ is the first invariant of the stress tensor, $c$ is the cohesion which may depend on the deformation history; $\alpha_1$ and $\alpha_2$ are material-dependent constants which are functions of the internal friction angle ($\varphi$) as follows for
the specific case of a plane strain deformation (de Souza Neto et al. 2011):

$$ \alpha_1 = \frac{3 \tan \varphi}{\sqrt{9 + 12 \tan^2 \varphi}}, \quad \alpha_2 = \frac{3}{\sqrt{9 + 12 \tan^2 \varphi}}. \qquad (12) $$

Viscoplastic flow rate initiates when rate-dependent plastic deformation has begun and follows the standard formulation (Desai & Zhang 1987; Vermeer 1990; Abu Al-Rub & Tehrani 2011):

$$ \dot{\varepsilon}_{ij}^{\text{vp}} = \dot{\gamma}^{\text{vp}}(\sigma_{ij}^{\text{v}}, T) \frac{\partial \Phi_{F}^{\text{vp}}}{\partial \sigma_{ij}^{\text{v}}} = \dot{\gamma}^{\text{vp}}(\sigma_{ij}^{\text{v}}, T) \frac{\partial}{\partial \sigma_{ij}^{\text{v}}} \left( \sqrt{J_{II}} + \alpha_3 P \right) \quad (13) $$

specifying a magnitude and direction of viscoplastic flow, with $\dot{\gamma}^{\text{vp}}$ as a non-negative quantity defined as the viscoplastic consistency parameter describing the magnitude of viscoplastic flow (Simo & Hughes 1998); $\Phi_F^{\text{vp}}$ is the viscoplastic flow potential and $\alpha_3 = 3 \tan \psi / \sqrt{9 + 12 \tan^2 \psi}$, $\psi$ is the dilatancy angle whose value may describe associative or non-associative plasticity (Vermeer & de Borst 1984). The derivative of the viscoplastic flow potential with respect to the stress tensor describes the direction of viscoplastic flow:

$$ \frac{\partial \Phi_F^{\text{vp}}}{\partial \sigma_{ij}^{\text{v}}} = \frac{s_{ij}^{\text{v}}}{2\sqrt{J_{II}(s_{ij}^{\text{v}})}} + \frac{\alpha_3}{3} \delta_{ij}. \qquad (14) $$

The equation above accounts for both deviatoric and dilatant viscoplastic flow. Non-zero dilatancy removes the assumption of a plastically incompressible material, that is, one for which pressure changes due to viscoplastic deformation do not result in a net volume change in the material (Poliakov & Herrmann 1994; Poliakov et al. 1994; Turcotte & Schubert 2002). In this crucial respect, our work differs from widespread practice in geodynamics, which approximates rocks as either elastically or plastically incompressible (Babeyko & Sobolev 2008; Schmeling et al. 2008; Kaus 2010; Leng & Gurnis 2015; Ruh et al. 2015; Jaquet et al. 2016; Kiss et al. 2020; Pajang et al. 2021).

The viscoplastic consistency parameter $\dot{\gamma}^{\text{vp}}$ can represent various functional forms depending on the problem and material (de Souza Neto et al. 2011), we utilize the functional form of (Zienkiewicz & Cormeau. 1974; Owen & Hinton 1980; Perzyna 1986; Desai & Zhang 1987):

$$ \dot{\gamma}^{\text{vp}}(\sigma_{ij}^{\text{v}}, T) = \frac{1}{\mu} \left\langle \frac{\Phi_{Y}^{\text{DP}}}{\Phi_{0}} \right\rangle^{m}, \qquad (15) $$

where $\langle \cdot \rangle$ is the Macaulay bracket defined for any function f as: $\langle f \rangle = f$ if $f > 0$ and $\langle f \rangle = 0$ if $f \le 0$. m here represents a material parameter following the form of Desai & Zhang (1987). $\Phi_Y^{\text{DP}}$ represents the Drucker–Prager yield criterion, while $\Phi_0$ is a normalizing term usually taken as the yield stress or cohesion (Desai & Zhang 1987), or taken as the plastic viscosity (Jacquey & Cacace 2020). Considering viscoplasticity, the model requires an additional term, $\mu$ in eq. (15), where $1/\mu$ has units of the inverse of time and expresses the relative rate of viscoplastic strain (Desai & Zhang 1987). The temperature dependence of plastic deformation is included in the yield function $\Phi_Y^{\text{DP}}$ through $s_{ij}^{\text{v}}$. In principle, we can include frictional hardening (or softening) and cohesion hardening (or softening) depending on the deformation history through a variety of strategies (Leroy & Ortiz 1989, 1990), but these will not be considered here for the sake of simplicity.

### 2.2.4 Sources of volumetric strains

Volumetric strains measure volume changes in a material (Bower 2009). They are often suppressed in geodynamic codes which assume incompressibility by setting the divergence of velocity to zero (Kaus 2010; Thielmann & Kaus 2012; Duretz et al. 2015; Jourdon et al. 2018; Bessat et al. 2020). Since the volumetric strains are related to the pressure through the bulk modulus, setting it to an arbitrarily high value can suppress volumetric elastic strains (Zienkiewicz & Taylor 2005). These are volume changes from elastic rheologies. In terms of plasticity, not all flow laws accommodate dilatant behaviour, for example, the von Mises criterion only utilizes the deviatoric stress invariant and cohesion to define the plastic strength of a material; other flow laws can suppress the influence of dilatant strain when the dilatancy angle is set to zero (Babeyko & Sobolev 2008), thereby suppressing volume changes due to plastic flow.

In addition to accounting for elastic compressibility in our constitutive laws, we include contributions to volumetric strains from dilatant plasticity during irreversible brittle deformation. Our point will be to show specifically the impact of this brittle plastic compressibility (dilatant plastic strain) on the deformation state and thermal feedback.

## 2.3 Numerical implementation

For computational purposes, we express eqs (4) and (5) in incremental form:

$$ \Delta \varepsilon_{ij} = \Delta \varepsilon_{ij}^{e} + \Delta \varepsilon_{ij}^{v} + \Delta \varepsilon_{ij}^{vp}, \qquad (16) $$

$$ \Delta \varepsilon_{ij} = \Delta \varepsilon_{ij}^{e} + \Delta \gamma^{v} \frac{\partial \Phi_{F}^{v}}{\partial s_{ij}} + \Delta \gamma^{vp} \frac{\partial \Phi_{F}^{vp}}{\partial \sigma_{ij}}. \qquad (17) $$

Similarly, the unknowns in eqs (9) and (15) which were expressed in rate form are expressed in incremental form.

In solving the conservation equations, we have utilized a solid mechanical finite element solver, Abaqus® (Abaqus 2010, 2019). Abaqus® is a robust, optimized, engineering industry-accepted finite element solver with a wide array of rheologies and element library. The software package also allows a user to custom any rheology or process through subroutines. Whenever stresses are updated through customized subroutines, a consistent algorithmic tangent modulus is required to aid convergence and estimate a time increment for the next simulation time step. Abaqus® has been widely used to study a variety of geodynamic problems (Regenauer-Lieb & Yuen 1998; Branlund et al. 2000, 2001; Regenauer-Lieb et al. 2001; Gerbault et al. 2002; Regenauer-Lieb & Yuen 2003; Dyksterhuis et al. 2005; Capitanio et al. 2007; Salomon 2018). We have used Abaqus® capabilities to solve the mechanical conservation equations, while we implemented the constitutive laws and the energy conservation problem through our own customized subroutines. We therefore transcribed our mechanical and thermal constitutive laws described in Section 2.2 into Fortran subroutines: a User-Material (UMAT) and User-Material Heat Transfer (UMATHT), respectively, implementable in Abaqus®.

The details of the algorithmic implementation of our constitutive laws, including the detailed computation of the consistent algorithmic tangent moduli for viscoelastic and viscoplastic rheologies, are given in Appendix A. A summary of the thermomechanical implementation is illustrated in a pseudo-code in Algorithm 1. Our implementation can be summarized as an elastic trial state, a thermally activated viscous correction to the elastic state (viscoelastic
rheology), followed by a frictional plasticity correction activated upon the material stress state above the Drucker–Prager yield criterion.

**Algorithm 1.** Pseudo-code for custom viscoelastic-viscoplastic constitutive law, over one time-step loop.

**Require:** $\Delta\varepsilon_n^{\text{e, trial}}$ at time step n

From $\Delta\varepsilon_n^{\text{e, trial}}$, obtain $\sigma_n^{\text{e, trial}}$

if $\sqrt{\varepsilon_{\text{II}}^{\text{e}}} \neq 0$ then

Start creep routine

while $\Phi^{\text{v}} > \text{Tolerance do}$

Newton-Raphson iterations to obtain $\Delta\gamma^{\text{v}}$

end while

Update $s_n^{\text{v}}$ and $J_{\text{II}}^{\text{v}}(s_n^{\text{v}})$

Compute $\Phi_{\text{v}}^{\text{DP}}(s_n^{\text{v}})$

if $\Phi_{\text{v}}^{\text{DP}} > 0$ then

Viscoplastic routine

while $\tilde{\Phi}_{\text{R}}^{\text{vp}} > \text{Tolerance do}$

Newton-Raphson iterations to obtain $\Delta\gamma^{\text{vp}}$

end while

Viscoplastic updates

Update $\varepsilon_{n+1}^{\text{vp}}$, $s_{n+1}^{\text{vp}}$, $P_{n+1}$, $\sigma_{n+1}$, $T_{n+1}$.

Consistent algorithmic tangent modulus.

Compute $C_{ijkl}^{\text{vp}}$

else

Viscous updates

Update $\varepsilon_{n+1}^{\text{v}}$, $s_{n+1}^{\text{v}}$, $\sigma_{n+1}$, $T_{n+1}$.

Consistent algorithmic tangent modulus

Compute $C_{ijkl}^{\text{v}}$

end if

else

The stress state is elastic

$$\sigma_{n+1} = \sigma_n^{\text{e, trial}}$$

end if

At each simulation time step for all integration points (element nodal), Abaqus calls our UMAT subroutine to compute the stresses and uses the stresses to solve the mechanics (eq. 1). The consistent algorithmic tangent modulus is then used to estimate the adaptive simulation time step to be used for the next time increment.

In coupling the mechanical deformations (eq. 1) to the energy conservation problem (eq. 3), we compute the heat source given by the second term on the right hand side of (eq. 3), in the UMAT after updating the stresses which is then passed to the User-Material Heat Transfer (UMATHT) Fortran subroutine. Finally, the stresses, temperature and other state variables are stored, which are called at the next time step as history-dependent variables.

Since we are interested in the feedback from mechanical work to thermal dissipation and vice versa, we utilized the 2-D coupled continuum plane strain temperature-displacement triangular element library (CPE3T) with linear shape functions within the Abaqus® element library. The CPE3T element library is composed of three nodal points, with two of them having displacement degrees of freedom and one temperature degree of freedom.

In addition to accounting for dilatant heating terms in the energy budget, another originality of our approach is that we utilize two rounds of Newton-Raphson iterations to compute the viscoelastic and viscoplastic multipliers. To ensure a strict positivity of these terms for the respective rheologies and convergence of the local Newton-Raphson iterations, we include a simple bisection scheme

<figure><img src="image_2.png" /><figcaption>**Figure 1.** TEST 1, input model for studying crustal-scale shear-banding in an isothermal medium with a rock matrix characterized by Maryland diabase rheology and a semicircular weak zone of radius 3 km characterized by Westerly granite. The setup was introduced by Duretz *et al.* (2014) to study the impact of viscous shear heating on the thickness of shear zones. Velocities are imposed on the top and lateral boundaries to satisfy a background strain rate of $5 \times 10^{-14} \text{ s}^{-1}$; all boundaries have zero shear stresses, and the bottom of the model is free-slip.</figcaption></figure>

within the Newton-Raphson loops (de Souza Neto 2004; Chapra & Canale 2011).

## 3 ALGORITHMIC TESTING AND VALIDATION

To assess the performance of our constitutive laws, we investigated shear band formation using viscoelastic rheology and then, viscoelastic-viscoplastic rheology. We first proceeded by investigating the formation of ductile shear bands using viscoelastic rheology as in Duretz *et al.* (2014) with an initial constant temperature. Thereafter, we investigated crustal-scale shear band formations in brittle and ductile regimes using viscoelastic-viscoplastic rheology as used by Duretz *et al.* (2021). We utilized the same geometry and boundary conditions for a systematic comparison with published results and highlighted similarities and differences in our different constitutive approaches.

### 3.1 TEST 1 (benchmark test 1): localization in an isothermal viscoelastic medium

#### 3.1.1 Model configuration and boundary conditions

Here, we utilized the setup introduced by Duretz *et al.* (2014). The model setup has dimensions of 70 km by 40 km, comprises a rock matrix that approximates a Maryland diabase rheology and a 3 km radius semicircular weak inclusion which approximates dry Westerly granite rheology centred on the bottom boundary of the model as shown in Fig. 1. The material properties for the matrix and weak inclusion are shown in Table 1.

The model boundary conditions approximate a pure shear experiment; therefore we apply velocities on the lateral and top boundaries to satisfy a constant background strain rate of $5 \times 10^{-14} \text{ s}^{-1}$, that is, a time-varying velocity boundary condition; while the bottom boundary is free to slip. All boundaries are thermally insulated so that no heat is evacuated to the surroundings or admitted therefrom. The model was discretized with 20 350 triangular elements connected by 10 477 nodes with linear shape functions, accounting for displacement and temperature degrees of freedom. The element type we use is coupled thermal and displacement elements in plane
<table><caption>Table 1. TEST 1 rheological parameters used in the benchmark setup shown in Fig. 1 drawn from Duretz et al. (2014).</caption><thead><tr><th>Material</th><th><i>A</i> (Pa<sup>-m</sup>s<sup>-1</sup>)</th><th><i>m</i></th><th><i>E</i><sub>a</sub> (kJ mol<sup>-1</sup>)</th><th>α<sub>th</sub> (m<sup>2</sup>s<sup>-1</sup>)</th><th><i>C</i><sub>p</sub> (J kg<sup>-1</sup>K<sup>-1</sup>)</th><th>ρ (kgm<sup>-3</sup>)</th></tr></thead><tbody><tr><td>Matrix</td><td>3.20×10<sup>-20</sup></td><td>3</td><td>276</td><td>8.82 × 10<sup>-7</sup></td><td>1050</td><td>2700</td></tr><tr><td>Inclusion</td><td>3.16×10<sup>-26</sup></td><td>3.3</td><td>186</td><td>8.82 × 10<sup>-7</sup></td><td>1050</td><td>2700</td></tr></tbody></table>

strain configuration which can help in the efficient conversion of mechanical work to heating provided the heat source is estimated.

As Duretz et al. (2014) used viscous rheology for an incompressible fluid with viscous feedback to heating (i.e. heat source term described by the inner product of the viscous deviatoric stresses and the viscous strain rate tensors), we utilized the viscoelastic approach with viscous thermal feedback for our test. We utilized a Young’s modulus of 25 GPa. We did not enforce a zero velocity divergence condition as is done for incompressible flow (Duretz et al. 2014). While the power-law rheology is formulated in terms of effective viscosity in Duretz et al. (2014), we implement our creep rheology in terms of strain rate. An additional difference is that we utilized a finite element approach allowing the elements to deform with the model, while Duretz et al. (2014) utilized the finite difference/marker in cell approach (Gerya & Yuen 2003). Finally, we adapted our time stepping based on the deformation state at the preceding time step, while fixed time steps were used in the reference study (Duretz et al. 2014).

### 3.1.2 Results of viscoelastic benchmarking

Localized deformation originated from the weak inclusion and propagated symmetrically towards the top left and top right corners of the model domain. In terms of the geometry of the model compared to the benchmark setup (Duretz et al. 2014), we were able to realize the same amount of bulk shortening (Fig. 2a). The logarithm of the second invariant of strain rate tensor, $\sqrt{\dot{\varepsilon}_{ij}\dot{\varepsilon}_{ij}/2}$, is shown in Fig. 2(a) after an estimated 25 per cent of bulk shortening. This evolution of strain rate reproduces the results of the reference model (Duretz et al. 2014) in terms of the amplitudes and width of the shear bands. We note some smearing as the shear bands propagated towards the boundaries. Shown in Fig. 2(b) is the temperature accumulated at the end of the simulation due to viscous shear heating. This pattern aligns with the viscoelastic deformation, which is the source of the dissipative heating. Temperature rose by up to 205 K within the shear band.

We extracted a 1-D plot of the difference between the final temperature and the initial temperature in an element within the shear band near the top left corner of the deformed model of Fig. 2(b), and compared the evolution with that of Duretz et al. (2014) as shown in Fig. 2(c). The temperature perturbation is similar to that previously obtained by Duretz et al. (2014) until about 15 per cent of axial shortening, after which a difference of 5 K builds up at 25 per cent of axial shortening. However, this difference of 5 K between our results and the reference study only represents about 3 per cent given the large temperature increase. We attribute this difference to differences in the numerical schemes mentioned previously, which is further supported by the supplementary tests carried in the following section.

### 3.1.3 Sensitivity to different element types, linear and nonlinear interpolators

In order to assess the mesh sensitivity of our viscoelastic implementation, especially the influence of mesh types and the type of

interpolation used, we used the set up shown in Fig. 1. For the results presented in Fig. 2(a), we used linear triangular elements in which the displacements between nodes were interpolated linearly. In assessing the robustness of our development, we carried out resolution tests using coarser (7468) and finer (42 604) triangular elements compared to the elements used in Fig. 2(a). The resolution tests indicate that the shear bands are resolved irrespective of the resolution (Figs 3a and b). Using different element types (quadrilateral elements) and including nonlinear (quadratic) shape functions, which interpolate displacements between nodes using higher order polynomials, we also found that the deformation was not mesh sensitive (Figs 3c and d).

We utilized the same setup as shown in Fig. 1 but with 20 411 quadrilateral elements, with linear and nonlinear shape functions. As shown in Fig. 3, the dimensions of the deformed domain as well as the width of the major shear bands are essentially identical to those obtained using triangular elements (Fig. 2a). The results are, however, insensitive to whether we used low-order or high-order finite elements. Apart from slightly higher strain rates when quadrilateral elements were used compared to triangular elements, the width of the shear bands and the associated temperature were insensitive to the different mesh types, interpolating functions and resolution.

## 3.2 TEST 2 (benchmark test 2): localization in brittle and ductile regimes

### 3.2.1 Model setup and boundary conditions

The setup discussed for viscous shear-banding in a ductile regime in Section 3.1 above corresponds to an isothermal regime without body forces and with viscoelastic rheology to conform with the reference model (Duretz et al. 2014). To test a case incorporating our viscoelastic-viscoplastic rheology, we utilize another setup introduced by Duretz et al. (2020, 2021); it assumes a 100 km long and 30 km thick crust with rheology corresponding to that of Westerly granite, and in which a 2 km radius circular imperfection is embedded in the middle of the model, as shown in Fig. 4. The model includes an initial temperature of 293 K at the top increasing linearly to 739 K at the base of the model (15 K km<sup>-1</sup>). All boundaries are insulated to prevent heat loss to the surroundings or input from external sources, such that dissipative heating only arises from the internal deformation. The material properties are shown in Table 2. We carried out simulations to keep a constant boundary strain rate of 10<sup>-15</sup> s<sup>-1</sup>. Note that we used an additional parameter μ described in eq. (15), where 1/μ is the relative rate of viscoplastic strain, chosen to be 10<sup>-15</sup> s<sup>-1</sup>.

Points of departure from the results of Duretz et al. (2021) and our model arise from several differences in the numerical methods, which were mentioned already above. Furthermore, in Duretz et al. (2018, 2021), the conservation equations were discretized using a staggered grid finite difference scheme where velocity and pressure are considered as primary variables, as traditionally assumed by the geodynamic community; but here, we have discretized our system of equations on a finite element grid with displacement as primary
<figure><img src="image_3.png" /><figcaption>**Figure 2.** TEST 1 results: (a) viscoelastic strain rate, (b) accumulated temperature shown in deformed configuration. The bulk shortening is 25 per cent and a background strain rate of $5 \times 10^{-14} \text{ s}^{-1}$. (c) Temperature increase for an element within the shear band location of element indicated by an arrow in Fig. 2(b) compared to temperature evolution within the shear zone of Duretz *et al.* (2014).</figcaption></figure>

variable, following a typical engineering approach (Zienkiewicz & Taylor 2005).

### 3.2.2 Strain localization in brittle and ductile regimes

The models were run for $1.2811 \text{ Myr} (4.04 \times 10^{13} \text{ s})$ with an adaptive time stepping. The initial boundary strain rate and the time of the experiment assure a shortening of only 4 per cent, consistent with our assumption of small strain theory. We carried out simulations with zero dilatancies without heating (Fig. 5a), 10° dilatancy without heating (Fig. 5b), with shear heating (Fig. 5c) and with shear and dilatant heating (Fig. 5d). Whether viscoelastic or viscoplastic deformation dominated at any time step, both regimes were characterized by the formation of shear bands initializing within the weak zone and propagating into the brittle and the ductile domains.

The upper brittle domain conformed to a viscoplastic deformation, characterized by narrow shear bands in all simulation cases;

while the lower ductile domain conformed to a viscous deformation, characterized by comparatively broader shear bands which propagated towards the bottom boundaries (Figs 5a–d). While the shear bands attenuated in both upper and lower domains towards the boundaries, the attenuation occurred faster in the hotter ductile regions. The transition from brittle to ductile behaviour was observed within a depth range of 14.21 and 17.7 km, characterized by a distinct horizon and peak pressure, except within the initial weak zone where pressure remains minimal (Figs 5a–d).

Our result with $\psi = 10^°$ succinctly captured similar features to the work of Duretz *et al.* (2021) with 10° dilatancy angle and bulk modulus of $5 \times 10^{10} \text{ Pa}$, despite using different numerical discretization schemes and constitutive law updates. These features include: (1) the vertical symmetry of the deformation in terms of the propagation and position of shear bands (upper versus lower domains), (2) the style of dominant shear bands (narrower in the brittle regime and broader in the ductile regime), (3) ~33° angle of shear bands compared to the range of ~30° and 35° of the reference
<figure><img src="image_4.png" /><figcaption>**Figure 3.** TEST 1, resolution tests and element type with linear and nonlinear interpolators: (a) logarithm of the second invariant of strain rate for a coarse and finer triangular mesh compared to Fig. 2(a), (b) corresponding temperature for the coarse and finer meshes, (c) logarithm of second invariant of strain rate invariant using linear (low order) and nonlinear (high-order) 20 411 quadrilateral elements and (d) corresponding temperature for using low-order and high-order elements (linear versus quadratic interpolants). This experiment was done for 25 per cent of shortening.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>**Figure 4.** TEST 2: input model configuration for studying crustal-scale shear-banding in a pure-shear compression experiment with a rock matrix characterized by Westerly Granite rheology and a circular weak zone of radius 2 km at the centre. Velocities are imposed as shown to satisfy a background of $10^{-15}$ s$^{-1}$. The temperature at the top is 293 K with a linear gradient to 739 K at the bottom. We include a pre-loading gravity stage before imposing the velocity boundary conditions. The top surface is kept as a free surface throughout. The model setup and boundary conditions were introduced by Duretz *et al.* (2020, 2021).</figcaption></figure>

<table><caption>**Table 2.** TEST 2: rheological parameters used in the benchmark setup shown in Fig. 4. Bulk modulus, $K = 5 \times 10^{10}$ Pa, $C_p = 1050$ J kg$^{-1}$K$^{-1}$, $\rho = 2700$ kgm$^{-3}$ and a background strain rate was set at $10^{-15}$ s$^{-1}$. The dilatancy angle $\psi$ (°) ranges between 0 and $10^\circ$ for the matrix and was kept at $0^\circ$ for the inclusion (see the text for details). These parameters were drawn from Duretz *et al.* (2020).</caption><thead><tr><th>Material</th><th>A (Pa<sup>-m</sup>s<sup>-1</sup>)</th><th>m</th><th>E<sub>a</sub> (kJ mol<sup>-1</sup>)</th><th>α<sub>th</sub> (m<sup>2</sup>s<sup>-1</sup>)</th><th>φ (°)</th><th>c<sub>0</sub> (MPa)</th></tr></thead><tbody><tr><td>Matrix</td><td>3.16×10<sup>-26</sup></td><td>3.3</td><td>186.5</td><td>8.82 × 10<sup>-7</sup></td><td>30</td><td>50</td></tr><tr><td>Inclusion</td><td>1×10<sup>-20</sup></td><td>1</td><td>0</td><td>8.82 × 10<sup>-7</sup></td><td>0</td><td>0.1</td></tr></tbody></table>

model (Duretz *et al.* 2021), corresponding to theoretical Arthur or intermediate mechanically stable shear band angles defined by $45^\circ \pm (\varphi + \psi)/2$ (Arthur *et al.* 1977; Kaus 2010) and (4) the brittle-ductile transition (just below the weak inclusion), amplitudes of deformation as well as the pressure evolution (Fig. 5c).

Including heat production (either shear heating or shear and dilatant heating were accounted for) smoothes out the shear bands in the brittle and ductile domains (Figs 5c and d). In terms of 1-D comparisons of volumetric strain rate (calculated as $\dot{\varepsilon}_{kk}$) and deviatoric strain rate invariant, we observe that they are not affected by either shear heating or shear and dilatant heating (Figs 5e and f), except that dilatant strain rates are higher for non-zero dilatancy and marginally higher for $10^\circ$ dilatancy without heating. Localization is higher for non-zero dilatancy (Fig. 5f). We also observed a temperature increase of between 31.42 and 35.32 K for shear heating and 29.73 and 32.66 K for shear and dilatant heating corresponding to a temperature reduction of between 1.69 and 2.66 K. The difference between the second invariant of stress between zero-dilatancy and $10^\circ$ dilatancy is between 16.38 and 39 MPa, while the pressure difference ranged between 30 and 75.2 MPa (Figs 5e and f). Despite the constant background strain rate in these simulations, dilatant plastic feedback through dilatancy has a small effect on the temperature evolution.

### 3.3 The role of plastic dilatancy

Since the contribution to volumetric stresses and strains includes dilatant plasticity, we investigated the effect of varying the dilatancy angles, for a given value of elastic compressibility ($K = 5 \times 10^{10}$ Pa). We set the activation energy to zero for the weak inclusion to ensure viscous behaviour therein irrespective of temperature (Table 2). The results for constant dilatancy angles $\psi = 0^\circ$ to $30^\circ$ are shown in Fig. 6 for volumetric and deviatoric strain rate invariants,

and Fig. 7 for deviatoric stress invariant and pressure.

While the volumetric strain rate invariant increased in amplitude from non-dilatant plasticity to the extreme case of associated plasticity with dilatancy angle $= friction angle = 30^\circ$, we observed uniform volumetric strains in the ductile domain (below 15 km from the surface) as shown in Fig. 6. Concerning the deviatoric strain rate invariant, the shear bands are narrow in the brittle domain (above 15 km) and broader in the ductile domain (below 15 km), as shown in Fig. 6. Increasing dilatancy angles leads to increasing strain magnitudes within the shear bands in the brittle domain.

The evolution of the deviatoric stress invariant and the pressure indicate a gradual increase in both quantities near the brittle-ductile transition as the dilatancy angles increased, with the pressure saturating around the brittle-ductile transition for dilatancy of $30^\circ$ (Fig. 7). The zone of highest deviatoric stress invariant and pressure also moved upwards as dilatancy angles increased.

To further assess the evolution of deformation and stress state, we compared the evolution of the volumetric strain rate invariant and deviatoric stress invariant, we extracted 1-D curves for both quantities (Fig. 8). Due to the vertical symmetry of deformation, we displayed plots at the location indicated in Fig. 6(a) and Fig. 7(a). Concerning the volumetric strain rate invariant, we observed that increasing dilatancy angles leads to an increase in the volumetric strain rate invariant in the brittle domain (Fig. 8a). With respect to the second invariant of the stress tensor ($J_{II}$), the curves intersect near the brittle-ductile transition below 15 km depth (Fig. 8b). A first-order observation is the similarity in the shapes of the curves. The maximum deviatoric stress invariant ranged from 531 MPa for $\psi = 0^\circ$ and 618 MPa for $\psi = 30^\circ$ at depths of 15.4 and 13.9 km, respectively (Fig. 8c). This indicates that higher dilatancy angles require higher stress at a relatively shallower depth. Even though the experiment is a pure shear case, the impact of dilatancy is more obvious in the brittle domain. Decreasing dilatancy angles from realistic values of $10^\circ$, that is, at least $20^\circ$ less than the friction angle (Vermeer & de
<figure><img src="image_6.png" /><figcaption>**Figure 5.** TEST 2 results shown in undeformed configuration from the benchmark setup of 100 km × 30 km (Fig. 4). The left panel from (a) to (d) shows the logarithm of the second invariant of strain rate; while the right panel indicates the pressure. Model with: (a) zero dilatancy without dissipative heating, (b) 10° dilatancy without dissipative heating, (c) 10° dilatancy with shear heating, (d) 10° dilatancy with shear and dilatant heating. The arrow at 7.5 km below the surface of the model in (a) indicates the starting location where horizontal profiles were extracted: (e) volumetric strain rate, (f) logarithm of deviatoric strain rate, that is, second invariant of strain rate, (g) temperature without heating, with shear heating alone and shear plus dilatant heating and (h) pressure and the second invariant of deviatoric stress.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Figure 6.** The effect of constant plastic dilatancies on volumetric strain rate invariant (left half panel from $x = 0$) and deviatoric strain rate invariant (right half panel from $x = 0$) where shear and dilatant plastic deformations contribute to dissipative heating. The dilatancy angles span (a) 0° (non-dilatant plasticity) to (g) 30° (an extreme case of the friction angle = dilatancy angle). The black arrow indicates the location where 1-D profiles were extracted and shown in Fig. 8. We have shown one half of each quantity due to the vertical symmetry about $x = 0$.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Figure 7.** The effect of constant plastic dilatancies on deviatoric stress invariant (left panel) and pressure (right panel) where shear and volumetric plastic deformations contribute to dissipative heating. The dilatancy angles span (a) 0° (non-dilatant plasticity) to (g) 30° is (an extreme case of the friction angle = dilatancy angle). The black arrow in (a) indicates the location where 1-D profiles were extracted and shown in Fig. 8.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Figure 8.** 1-D plots of (a) volumetric strain rate invariant and (b) deviatoric stress invariant with locations indicated by arrows in Figs 6(a) and 7(a), (c) plot of the maximum deviatoric stress invariant (used as the limit of the brittle domain) and depth to the brittle-ductile transition with respect to the dilatancy angles.</figcaption></figure>

Borst 1984) to non-dilatant plasticity led to a 25 MPa (or 10 per cent) reduction in the deviatoric stress invariant near the surface, and ~30 MPa (or 5 per cent) reduction at ~14.5 km before the brittle-ductile transition.

If we use this second stress invariant as a proxy for the magnitude of shearing exerted by a state of stress in the medium (Bower 2009), we may argue that lower dilatancy angles reduce the magnitude of shear stresses in the medium compared to higher dilatancy angles. The synchrony of the deviatoric stress invariants for different dilatancy angles in the ductile domain can be attributed to the assumption that there are no volumetric plastic deformations in the ductile domain, per our constitutive description. The apparent attenuation in the second stress invariant for small dilatancy angles may explain an increase in strain rate, alluding to the inverse relationship between the strain rate and second stress invariant (eq. A9).

## 4 THE EFFECT OF VOLUMETRIC PLASTIC (DILATANT) HEATING

The pure shear experiment was implemented with the deliberate choice of constant strain rate boundary conditions by adapting the boundary velocity at each time step, and reproducing the conditions for benchmarking purposes (Duretz et al. 2021). We showed that the contribution of dilatant plastic dissipation to volumetric or deviatoric localization was not significant (Figs 5e and f) and impacted the thermal state by reducing the temperature by up to 3 K (Fig. 5g).

To investigate how different boundary conditions may influence deformation and heat production when dilatant plastic dissipation is included, we ran another series of simulations: where we imposed constant boundary velocities (TEST 3), utilized an isothermal domain (TEST 4A) and high-strain rate experiments (TEST 4B).
<figure><img src="image_10.png" /><figcaption>**Figure 9.** TEST 3: Results of deformation with dissipative feedback from shear heating and shear with dilatant heating. (a) volumetric strain rate invariant, (b) logarithm of the deviatoric strain rate invariant, (c) evolved temperature at the end of the experiment, (d) temperature evolution with time for an element within the shear band, (e) temperature evolution with time for an element within the weak inclusion, (f) vertical temperature profile near the weak inclusion crossing one of the shear band and (g) horizontal temperature profile in the brittle domain along the model. The red arrowheads in (a) are locations where the temperature histories were tracked within the shear band and weak inclusion, respectively; while the arrows indicate starting locations where 1-D temperature profiles were extracted. The left panel in Figs 9(a)–(c) represents results from shear heating, while the corresponding right panels show the results when dilatant heating is included.</figcaption></figure>

## 4.1 TEST 3: Nucleation zone in brittle-ductile domain

In TEST 3, we utilized the same input setup shown in Fig. 4 with the same material parameters. In addition, we replace the varying velocity boundary conditions with constant boundary velocities. We ran the simulations for 2.5622 Myr (8.08 × 10¹³ s) for a target shortening of 8 per cent. The results with shear heating alone, and shear with dilatant plastic feedback are shown in Fig. 9.

Shown in Figs 9(a) and (b) are the volumetric strain rate invariant and deviatoric strain rate invariant, respectively, for cases of shear heating alone and shear with dilatant heating. The corresponding temperatures at the end of the simulations are shown in Fig. 9(c). We again note that the deformation is symmetric (vertical symmetry about the initial weak zone depth) similar to the pure shear experiment reported in TEST 2, with similar temperature evolution in 2-D. While the volumetric strain rate invariant followed the shear bands in the brittle domain, we note the uniformity of the volumetric strain rate invariant in the ductile domain (Fig. 9a), which is consistent with our assumptions of treating volumetric plastic deformations only from brittle plasticity and not ductile creep. The observation

of shear bands in the brittle and ductile regimes for the deviatoric strain rate invariant is also consistent with our accounting for deviatoric strains in both the brittle and ductile domains. To compare the temperature, we extracted the temperature evolution for an element within a shear band and within the weak zone (Figs 9d and e). Within the shear band element, the temperature increased by 102 K for shear heating. We observed that both the shear heating and shear with dilatant heating temperatures evolved similarly until ~0.4 Myr where the temperatures were offset until the end of the experiment with the inclusion of dilatant heating leading to reduction of 5 K at the end of the simulation. In the case of the element within the weak inclusion, the overall temperature rise was 20 K at the end of the experiment when shear heating was considered alone; however, this temperature rise was reduced by 5 K when dilatant heating was included. The vertical temperature profile also shows that dilatant heating reduced the amount of heat dissipated by shear heating by 4 K (Fig. 9f). The same observation is true for the horizontal profile where the inclusion of dilatant heating reduced the heat produced by shear heating by 5.5 to 6 K (Fig. 9g). The temperature due to shear dissipative heating versus temperature due to shear and dilatant dissipative feedback indicates that
dilatant contributions introduce a less heating effect in the brittle domain.

## 4.2 TEST 4A: Viscoplasticity in an isothermal domain and possible rheological change due to thermal dissipation

Motivated by the impact of dilatant plastic deformation on thermal dissipation, we returned to the initial setup for viscoelastic benchmarking shown in Fig. 1. We utilized the material properties discussed for TEST 3, and we have included dilatancy angle of 10°. The initial temperature we have used was set at 473 K, to ensure that the impact of ductile creep was suppressed in order to ensure a viscoplastic initial rheology. We have run the simulations for $4 \times 10^{12}$ s (126 kyr).

We observed differences in the volumetric strain rate invariant (Fig. 10a), deviatoric strain rate invariant (Fig. 10b), the deviatoric stress invariant (Fig. 10c) and the temperature at the end of the simulations (Fig. 10d). The shear bands were observed at an angle of ~35–42°. Utilizing a friction angle of 30° and dilatancy of 10°, the angle of the observed shear bands within the matrix corresponds to Arthur or intermediate angles, suggesting that plasticity parameters influenced the behaviour of shear bands during viscoplastic flow as opposed to the ductile shear bands which are not influenced by these parameters.

When dilatant heating was included, shear bands were smoother as seen in the volumetric strain rate invariant and deviatoric strain rate invariant (Figs 10a and b). While the lowest deviatoric stress invariant was observed in the weak inclusion whether shear heating or shear and dilatant heating were considered, we observed a higher deviatoric stress invariant when dilatant heating was included (Fig. 10c). The temperature at the end of the simulation showed that temperatures increased for both shear heating and shear with dilatant heating (Fig. 10d). However, temperatures were higher within the shear bands when shear heating was considered alone compared to when dilatant heating was included (Fig. 10d). Our deduction from these results is that including the dilatant term contributed to a higher deviatoric stress invariant, and lower deviatoric strain rate invariant with less heating effect within the shear bands.

To further illustrate how dilatant plastic feedback may affect the behaviour of the plotted variables in Fig. 10, we extracted 1-D plots and history plots from locations shown in Fig. 10(a). Across the model, the contribution of dilatant plastic heating is either a reduction or an increase in the volumetric strain rate invariant (Fig. 11a) and the deviatoric strain rate invariant (Fig. 11b). We further confirmed that including dilatant heating led to increased pressure and deviatoric stress invariant in the brittle domain, with the pressure reducing within the shear bands for both shear heating and shear with dilatant heating between $x = -7.5$ and $-20$ km (Fig. 11c). This observation is consistent with our correction for pressure during viscoplastic flow (eq. A11). That is, zones with the highest strains (volumetric or deviatoric) follow from zones with the highest viscoplastic multiplier ($\Delta\gamma^{\text{VP}}$), hence the more pressure correction. The pressure difference ranges between 202 and 509 MPa when dilatant heating is included, while the difference in deviatoric stress invariant ranges is between 187.82 and 435 MPa. Across the profile, the temperature rise was between 28.9 and 126.6 K when shear heating was considered alone and between 24.5 and 99.6 K when dilatant heating was included with shear heating (Fig. 11d). The temperature reduction when dilatant heating was included was between 4.3 and 38 K (Fig. 11e).

In Figs 11(f)–(j), we show the history-dependent evolution of some of the variables for an element within the shear band whose location was shown in Fig. 10(a), near the weak zone. We observed that the integrated volumetric and deviatoric viscoplastic strain invariants indicate that dilatant heating combined with shear heating influences both quantities (Figs 11f and g). This indicates that at longer simulation times, dilatant heating may increase irreversible strains. This is consistent with the along-axis profile where it was observed that when dilatant heating was included, it led to an increase or decrease of volumetric strain rate or deviatoric strain rate as shown in Figs 11(a) and (b). The evolution of pressure and deviatoric stress invariant approached a steady state behaviour at ~83.6 kyr (Fig. 11h) at which point the inclusion of dilatant heating offset the temperature by about 18 K (Figs 11i–j). While the irreversible strains, pressure and stress invariants were similar until ~76 kyr, the temperature for either: shear heating alone and shear with dilatant heating were being offset from 20 kyr and the offset increased as the simulations proceeded.

It can be argued that shear heating alone can enhance deviatoric strain localization as already shown in previous studies (Thielmann & Kaus 2012), while shear with dilatant heating also enhances deviatoric strain localization even if delayed (Fig. 10b), increases deviatoric stress invariant (Fig. 10c) consistent with an inverse relationship between strain rate and deviatoric stress invariant (eq. A9) and reduce the amount of heat produced (Fig. 10d). Note that the same boundary conditions used here ensure a constant strain rate, which suggests that even for nominally volume-preserving boundary conditions, volumetric plastic effects contribute to the rate of deformation and temperature evolution depending on the rheology. It must also be noted that as temperature increases due to heating feedback, the likelihood of the rheology switching from viscoplastic to viscoelastic is not ruled out. In fact, in some parts of the model, this rheological switch was already seen when the material stopped yielding, that is, in these areas where the Drucker–Prager criterion (eq. 11) was no longer satisfied.

## 4.3 TEST 4B: Viscoelasticity and viscoplasticity at high strain rates

To perform experiments at high strain rates, we returned to the setup shown in Fig. 4 where we increased the initial boundary strain rate by increasing the initial velocity by three orders of magnitude to $1 \times 10^{-8} \text{ m s}^{-1}$ for the vertical boundaries and to $3 \times 10^{-9} \text{ m s}^{-1}$ for the bottom boundary compared to TEST 3. The total amount of shortening was 24 per cent. We ran the simulations for $1 \times 10^{12}$ s (~32 kyr) and show the results in Figs 12(a)–(c) and associated 1-D profiles for temperature variation across the brittle domain and temperature history within an element in the shear band are shown in Figs 12(d)–(g). The shear bands initiated within the weak zone and propagated through the brittle and ductile domain, secondary shear bands appeared on the top boundary as shown in Figs 12(a) and (b). We also observed that the ductile shear bands had higher amplitudes in the deviatoric strain rate invariant than brittle shear bands (Fig. 12b). We also followed the deformed zones where the heat is produced as shown in the temperature evolution (Fig. 12c) where more heat was produced in the ductile domain than the brittle domain.

In assessing the discrepancy in the deformation and the temperature when shear heating is considered alone or dilatant plastic contribution is incorporated, we observed from the horizontal profiles crossing two shear bands in the brittle domain that the feedback
<figure><img src="image_11.png" /><figcaption>**Figure 10.** TEST 4A results shown in deformed configuration for a setup shown in Fig. 1. We have utilized an initial isothermal temperature of 473 K to ensure viscoplastic deformations. The results show two columns indicating the influence of shear heating, and shear with dilatant heating, respectively: (a) volumetric strain rate invariant, (b) deviatoric strain rate invariant, (c) second invariant of stress and (d) temperature evolution at the end of the simulation. These results are plotted in deformed configuration measuring 57.4 km in the x-direction and 48.8 km in the y-direction similar to the same observation in Fig. 2. The arrows on the left border and one element within a shear band in (a) are pointers to, locations where 1-D profiles and some history-dependent variables were extracted and shown in Figs 11(a)–(e), and Figs 11(f)–(j), respectively.</figcaption></figure>

from dilatant and shear heating led to increased volumetric strain rate (Fig. 12d) and a slight increase or decrease in the deviatoric strain rate (Fig. 12e). Concerning temperature, the difference is up to 6 per cent (Fig. 12f), while an example of temperature evolution in an element within the shear band in the brittle domain shows that the dilatant plastic feedback can lead to an 8 per cent reduction in temperature after 30 000 yr. These results demonstrate, as with previous test cases, the ability of dilatant plastic dissipation to potentially slow down shear heating and impacts on the volumetric

strain rate. Over long time scales, this effect is expected to be more dominant.

## 5 APPLICATION TO LITHOSPHERIC DYNAMICS

To further illustrate the potential use of our proposed viscoelastic–viscoplastic constitutive law, we carried out simulations that tackle the lithospheric-scale.
<figure><img src="image_12.png" /><figcaption>**Figure 11.** TEST 4A (a–e) are 1-D profiles extracted across two dominant shear bands from location indicated by arrow on the left border of the model shown in Fig. 10(a): showing (a) volumetric strain rate invariant, (b) deviatoric strain rate invariant, (c) second invariant of stress shown in dashed lines and pressure shown in full lines, (d) temperature difference between the final temperature and initial temperature, and (e) final temperature; (f–j) are 1-D history evolution of quantities extracted from an element within one of the dominant shear bands whose location is shown in Fig. 10(a): (f) volumetric viscoplastic strain invariant, (g) deviatoric viscoplastic strain invariant, (h) deviatoric stress invariant shown in dashed lines and pressure shown in full lines, (i) evolution of temperature and (j) integrated temperature.</figcaption></figure>

## 5.1 Initial model set-up and boundary conditions

We investigated two models, both composed of a 12 km oceanic crust and a mantle: (1) without a weak zone, and (2) with an elliptical weak zone with major and minor radii 20 and 8.4 km, respectively, placed at 32 km below the surface. An example of the model set up with an elliptical weak zone is shown in Fig. 13. The model setup without the weak zone was aimed at investigating possible lithospheric buckling (or warping) when the model is subjected to long-term compression, while the model with the elliptical weak zone was to preferentially seed localization and to investigate the role of dilatant plastic and shear dissipation compared to shear dissipation alone. The material properties utilized are shown in Table 3.

The initial vertical thermal distribution in the mantle is computed using the cooling of a semi-infinite half-space model for a 50 Myr old lithosphere (Turcotte & Schubert 2002; Gerya 2019):

$$T_y = T_1 + (T_0 - T_1) \left( 1 - \frac{2}{\pi} \arctan \left[ \frac{y}{\sqrt{\alpha_{th} \tau}} \left( 1 + \left( \frac{y}{2\sqrt{\alpha_{th} \tau}} \right)^4 \right) \right] \right). \quad (18)$$

where $T_y$ is the temperature at a given depth, $y$, $T_0$ is the temperature at the surface, $T_1$ is the temperature at the base of the model; the function in parenthesis approximates the error function, where $\alpha_{th}$ represents thermal diffusivity and $\tau$ represents the age of the lithosphere in seconds. The initial temperature at the top of the model is set at 273 K, while the temperature at the base is set at 1673 K: Both
<figure><img src="image_13.png" /><figcaption>**Figure 12.** TEST 4B comparison of deformation states when accounting for shear heating alone and when dilatant plastic contributions were included (see Section 4.3). (a) volumetric strain rate invariant, (b) deviatoric strain rate invariant and (c) accumulated temperature. 1-D profiles extracted from ~7.5 km (indicated by an arrow on the border of the model shown in Fig. 12a) for: (d) volumetric strain rate invariant, (e) deviatoric strain rate invariant and (f) corresponding temperature across the profile. (g) temperature history for an element within the shear band in the brittle part of the model domain (location of the element tracking the temperature within a shear band is indicated by an arrow in Fig. 12a).</figcaption></figure>

the left and right boundaries have zero heat flux. The top boundary is a free surface, and the bottom boundary is maintained free slip. We apply normal inward constant velocities on both vertical boundaries of the model at a rate of 2 cm yr⁻¹ for at least 5 Myr (1.5768 × 10¹⁴ s). We meshed the uppermost part of the model with 2.5 km element size and increased the element size towards the base of the model.

## 5.2 Potential contribution of volumetric plastic deformation to lithospheric localization, and heating

### 5.2.1 *The buckling case*

The results of the experiment without the weak zone after 5.7 Myr of horizontal compression are shown in Fig. 14. Seven shear bands developed whether shear or shear and dilatant heating were accounted
<figure><img src="image_14.png" /><figcaption>**Figure 13.** Illustrative model setup and boundary conditions for two cases: case 1 without an initial weak zone, and case 2 with an elliptical weak zone. Only case 2 is shown in this sketch. The horizontal compression velocity is 2 cm yr⁻¹.</figcaption></figure>

<table><caption>**Table 3.** Material properties in the geodynamic model. *E* (Young's modulus), *v* (Poisson's ratio), *K* (bulk modulus), *G* (shear modulus), *ρ* (density), *g* (gravitational constant), *C*<sub>p</sub> (specific heat capacity), *α*<sub>th</sub> (thermal diffusivity), *A* (pre-exponential multiplier), *E*<sub>a</sub> (creep activation energy), *R* (molecular gas constant), *φ* (friction angle), *ψ* (dilatancy angle), *c*<sub>0</sub> (initial cohesion), *m* (stress exponent) and 1/*μ* (relative rate of viscoplastic strain). The references for material parameters were drawn from <sup>a</sup>Ranalli (1995), <sup>b</sup>Goetze & Evans (1979), <sup>c</sup>Hirth & Kohlstedt (2003) and Thielmann & Kaus (2012).</caption><thead><tr><th>Property (unit):</th><th>Mantle<sup>a,b,c</sup></th><th>Crust<sup>a</sup></th><th>Weak zone<sup>a,b,c</sup></th></tr></thead><tbody><tr><td><strong>Mechanical</strong></td><td></td><td></td><td></td></tr><tr><td><em>E</em> (GPa)</td><td>25</td><td>25</td><td>25</td></tr><tr><td><em>v</em></td><td>0.25</td><td>0.25</td><td>0.25</td></tr><tr><td><em>G</em> (GPa)</td><td>10</td><td>10</td><td>10</td></tr><tr><td><em>K</em> (GPa)</td><td>16.7</td><td>16.7</td><td>16.7</td></tr><tr><td><em>ρ</em> (kg.m<sup>-3</sup>)</td><td>3300</td><td>2950</td><td>3100</td></tr><tr><td><em>g</em> (m.s<sup>-2</sup>)</td><td>9.8</td><td>9.8</td><td>9.8</td></tr><tr><td><em>m</em></td><td>3.5</td><td>2.3</td><td>4</td></tr><tr><td><strong>Thermal</strong></td><td></td><td></td><td></td></tr><tr><td><em>C</em><sub>p</sub> (J kg<sup>-1</sup>K<sup>-1</sup>)</td><td>1050</td><td>1050</td><td>1050</td></tr><tr><td><em>α</em><sub>th</sub> (m<sup>2</sup>s<sup>-1</sup>)</td><td>8.7 × 10<sup>-7</sup></td><td>9.7 × 10<sup>-7</sup></td><td>9.2 × 10<sup>-7</sup></td></tr><tr><td><strong>Dislocation creep</strong></td><td></td><td></td><td></td></tr><tr><td><em>A</em> (Pa<sup>-1</sup>s<sup>-1</sup>)</td><td>2 × 10<sup>-21</sup></td><td>1.59 × 10<sup>-14</sup></td><td>2.08 × 10<sup>-23</sup></td></tr><tr><td><em>E</em><sub>a</sub> (kJ mol<sup>-1</sup>)</td><td>535</td><td>154</td><td>283</td></tr><tr><td><em>R</em> (J K<sup>-1</sup>mol<sup>-1</sup>)</td><td>8.31</td><td>8.31</td><td>8.31</td></tr><tr><td><strong>Plasticity</strong></td><td></td><td></td><td></td></tr><tr><td><em>φ</em> (°)</td><td>30</td><td>30</td><td>15</td></tr><tr><td><em>ψ</em> (°)</td><td>10</td><td>10</td><td>5</td></tr><tr><td><em>c</em><sub>0</sub> (MPa)</td><td>1</td><td>1</td><td>1</td></tr><tr><td>1/<em>μ</em> (s<sup>-1</sup>)</td><td>10<sup>-15</sup></td><td>10<sup>-15</sup></td><td>10<sup>-15</sup></td></tr></tbody></table>

for (Figs 14a and b). The interaction of these shear bands led to lithospheric buckling (Fig. 14c) and a zone of localization near the centre of the model (between *x* = 0 and −100 km). We also observed that volumetric and shear localization were more enhanced at the centre of the model when shear heating was considered alone, that is, the tent-shaped feature resulting from the intersection of the two central shear bands was stronger when shear heating was the only dissipation mechanism (Figs 14a and b). The temperature showed an uplift of isotherms where the shear bands interacted (Fig. 14d).

We compared the impact of shear dissipation and shear with volumetric plastic dissipation on the results of the buckling experiment. Fig. 15(a) displays the surface elevation for both cases with shear heating and with shear and dilatant heating. We extracted 1-D plots every 50 km along the model to investigate the variations of the magnitude of the vertical heat flow (Fig. 15b), depth to 895 K isotherm (Fig. 15c), depth to the 1473 K isotherm (Fig. 15d) and the depth to the maximum value of the second invariant of the stress tensor where the model switches from brittle behaviour to ductile

behaviour (Fig. 15e). The heat flow was calculated based on the temperature within the first 10 km using the standard formulation of the product of thermal conductivity and thermal gradient (England & Katz 2010). From the heat flow plots, we observed variations of up to 3 mW s⁻² reduction in heat flow when dilatant heating was included (Fig. 15b). We chose the 895 K isotherm as it was observed to be the isotherm that was elevated when shear heating was the dissipation mechanism (Fig. 14d). This isotherm was found to be offset by up to 2.5 km when volumetric plastic dissipation was included with shear dissipation (Fig. 15c). The depth to the 1473 K (ensuring viscous flow) indicates that volumetric plastic dissipation influences this isotherm (Fig. 15d). This suggests that including volumetric plastic dissipation thickens the brittle part of lithosphere by a few kilometres. Finally, the depth of the brittle-ductile transition is also influenced when the dilatant heating is included with shear heating (Fig. 15e). Considering that the maximum of the deviatoric stress invariant was observed to correspond to the depth of transition from brittle to ductile behaviour (Fig. 8b),
<figure><img src="image_15.png" /><figcaption>**Figure 14.** Geodynamic model results for setup in Fig. 13 without the weak zone after 5.7 Myr of horizontal compression corresponding to an overall shortening of 23 per cent. (a) volumetric strain rate invariant, (b) deviatoric strain rate invariant, (c) second invariant of stress and (d) evolved temperature.</figcaption></figure>

the results from the buckling experiment shown in Fig. 15(e) indicate that the depth to the brittle–ductile transition can be thickened by an additional 3.75 km near x = 0 when dilatant heating is included.

### 5.2.2 The domal case

Our experimental setup in Fig. 13 includes an elliptical weak zone which could represent a low-velocity zone associated with partial melting, internal fluid-filled fractures, high-temperature area or inherited features from a previous geological event. The importance of such weak features in influencing tectonic regimes has

been investigated by past authors (Cloetingh et al. 1982; Gurmis 1992; Nikolaeva et al. 2010; Ruh et al. 2015). Our interest here is to investigate if we can observe similar features associated with domal structures due to compressional tectonics, for example, an anomalous thermal zone beneath a dome, the formation of faults that may affect the surrounding areas and the surface topography.

For a simulation of 5.4 Myr, we show the deformed domain in an area of size 268.2 km by 600 km around the weak inclusion and away from the inevitable edge effects. Furthermore, since the structural evolution between cases accounting for shear heating and shear with dilatant heating were similar, we only show the results of
<figure><img src="image_16.png" /><figcaption>(b) Vertical heat flow</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>(c) Depth to 895 K isotherm</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>(d) Depth to 1473 K isotherm</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>(e) Depth to the maximum $\sqrt{J_{II}}$</figcaption></figure>

<figure><img src="image_20.png" /></figure>

**Figure 15.** Comparative profiles for geodynamic results shown in Fig. 14 (case 1) showing: (a) surface elevation, (b) magnitude of vertical heat flow, (c) depth to 895 K isotherm (d) depth to 1473 K isotherm and (e) depth to the maximum deviatoric stress where the model shifts from brittle to ductile behaviour.
<figure><img src="image_21.png" /><figcaption>**Figure 16.** Results for input model shown in Fig. 13 for shear and dilatant heating (case 2). (a) Logarithm of the deviatoric strain rate invariant, (b) temperature, (c) magnitude of vertical heat flow, (d) depth to 895 K isotherm chosen based on one of the isotherms that was shallower when dilatant heating was included; (e) depth to 1473 K where ductile behaviour was active, and (f) depth of the maximum deviatoric stress invariant which can be used as the transition to brittle mantle behaviour to ductile mantle behaviour.</figcaption></figure>

shear with dilatant heating for deformation (Fig. 16a) and temperature (Fig. 16b), where four shear bands emanated from the elliptical inclusion. Two of the shear bands propagated upwards to the brittle lithosphere and led to an elevated surface topography, while the other two propagated downwards into the ductile part before fading out (Fig. 16a). This behaviour led to the shallowing and deepening of some isotherms. We however compare results from shear heating and shear with dilatant heating in extracted profile plots which we show in figures (Figs 16c–f). The large-scale shear bands and the associated thermal feedback led to a zone with initially horizontal isotherms rising beneath the weak zone. While the dome is associated with elevated surface heat flow for shear or shear and dilatant heating (Fig. 16c), the depth to the 895 K isotherm exhibited a shallowing beneath the dome (Fig. 16d) and a deepening of the 1473 K isotherm (Fig. 16e). This suggests a thermal anomaly which is a consequence of shear and dilatant dissipation. The specific effect of dilatant heating is to shallow the 895 K isotherm by about 2 km beneath the dome. The depth to the maximum of the deviatoric stress invariant used as the point of transition from brittle mantle behaviour to ductile mantle behaviour is unaffected by shear heating or shear with dilatant heating (Fig. 16f). Comparing the results of the depths to the chosen isotherms and the depth to the

brittle–ductile transition, we can conclude that the effect of shear heating and shear with dilatant heating mostly influences the brittle domain in this case.

## 6 DISCUSSION

In this contribution, we have described a formulation for viscoelastic–viscoplastic (VE-VP) constitutive behaviour, incorporating elastic compressibility, plastic compressibility and power-law viscoplasticity. We addressed the potential feedback from volumetric plastic dissipation, which we referred to as dilatant heating, and had been traditionally neglected in conventional geodynamic modelling. Our approach therefore treated the domain as a solid–mechanical system which also differs from traditional approaches that utilize a fluid–mechanical idealization of a solid system. We have given a step-by-step description of the constitutive update for stresses and strains; consistent Jacobian matrices for adaptive time-stepping as well as a concise algorithmic description. In this concise discussion section, we shall examine the impact of this volumetric plastic feedback to the dissipative heating loop.
The justification for the effect of volumetric plastic deformation on heat reduction can be first explained in the following way: volumetric plastic deformation implies a work term whose sign depends ultimately on the constitutive laws. When thermal conductivity is low as is the case for rocks, the overall process should be close to isentropic, which implies that dilatant plasticity (expansion) indeed implies some cooling. From eq. (A11), once the viscoplastic multiplier is determined, a correction to the pressure is implemented. Since the pressure multiplier is always positive, there is a reduction in pressure.

For pure shear experiments (TEST 2), the impact of volumetric plastic and shear dissipation compared to shear dissipation alone did not show a large difference in the volumetric and deviatoric strain rates; pressure and deviatoric stress invariant. Despite the volume-preserving pure shear experiment, we observed a temperature reduction of up to 3 K when dilatant plastic heating was included. This is because shearing dominated the deformation, and the feedback to volumetric plastic deformations and resultant dilatant plastic heating was minimized, that is, there was little volume change, despite extensive internal shear deformations. This was indeed apparent in Figs 5(c) and (d), where the different dissipative feedback did not significantly impact the pressure.

TEST 4A to investigate viscoplasticity, which included an initially isothermal medium, indicated that volumetric plastic dissipation impacted the volumetric strain rate, deviatoric strain rate, pressure, deviatoric stress invariant and irreversible strains. History-dependent variables like the integrated irreversible strains were also influenced by the contribution of volumetric plastic dissipation. In all cases, the temperature rose with the feedback on the temperature commencing before the other history-dependent variables like accumulated viscoplastic strains and deviatoric stress invariants. High strain rate experiments (TEST 4B) also showed that temperature can be reduced by up to 25 to 50 K when dilatant plastic heating was considered.

It has been argued that rocks generally deform viscoelastically at low stresses and do not necessarily support elastic stresses over the timescale of a subduction process, for example (Schmeling et al. 2008). A metric for assessing the timescale over which elastic stresses can be relaxed is the Maxwell relaxation time, given as the ratio of the effective viscosity and the shear modulus. In some numerical studies, the impact of elasticity is inhibited by setting an unrealistically high value for the shear modulus, which makes the Maxwell relaxation time small, hence elastic stresses are dissipated quickly (Schmeling et al. 2008). We have not included assumptions in our model to mitigate elasticity, as it had been shown in other studies that elasticity (represented by the shear modulus) is an important attribute in lithospheric-scale strain localization (Jaquet et al. 2016; Bessat et al. 2020). Our inclusion of elasticity in our model is justified on the basis that, we assume that the material can deform elastically until thermally activated viscous flow kicks in, while viscoplastic behaviour evolves once the Drucker–Prager yield criterion is breached.

In utilizing our new implementation in studying lithospheric dynamics, our results indicated that realistic heat flow values were reproducible for thermal dissipation alone when compared with published results of heat flow measurements from compressional environments (Manga et al. 2012; Ezenwaka 2023). We also observed that the contribution of volumetric plastic dissipation either for the buckling experiment or for the collisional case was a reduction of heat flow values, and variation in the depth of some isotherms, with shear heating more likely to lead to lithospheric thinning than shear with dilatant plastic heating.

As observed in all the simulations reported here, it is much more common for several shear bands to develop (see Figs 14a and b, for example). This sets up the question of what happens where different shear bands intersect/interact. In these regions, deformation is both intense and going in all directions which impedes the development of a single dominant shear direction. Hence localization in the sense of allowing a simple macroscopic shear direction to develop is hampered. Nevertheless, deformational heating remains intense because strong irreversible deformation is still occurring in such regions (Fig. 14d).

## 7 CONCLUSION

We have developed a solid-thermomechanical constitutive description to study strain localization in compressional experiments with contributions from volumetric plastic dissipative feedback, which are usually not considered in geodynamic models. Using small strain assumptions, we developed and applied a full 2-D viscoelastic–viscoplastic stress update scheme and included time adaptivity with consistent algorithmic tangent moduli. Our constitutive laws have been implemented on the Abaqus® finite element solver. Based on our studies, we conclude as follows:

(i) We benchmarked our new solid-mechanical model against other published numerical models for viscoelastic and viscoelastic–viscoplastic rheologies. Our benchmarking tests indicate that our simulations are not mesh-sensitive nor are they sensitive to low-order (linear) or high-order (quadratic) elements.

(ii) In comparing results of shear heating and shear with dilatant heating for pure shear numerical experiments, our results indicated that the implicit conservation of volume in the experimental setup through the specific choice of boundary conditions to keep the background strain rate constant leads to the volumetric plastic feedback not having a strong influence on the results.

(iii) A departure from volume-conserving boundary conditions pure shear experiments demonstrated the contributions of volumetric plastic feedback to the temperature and the stress state. Using sequentially a linear temperature profile (TEST 3), an isothermal case that ensured viscoplastic deformation (TEST 4A) and a high boundary strain rate (TEST 4B), we observed that volumetric plastic dissipation contributed to the thermal state of our models by reducing the amount of heat produced, which affected the depth to isotherms and even the brittle-ductile transition.

(iv) Whether thermal dissipation is due to shear heating or shear with dilatant heating, both scenarios led to a temperature rise, but the inclusion of volumetric plastic dissipation led to reduced heating, sometimes up to 8 per cent.

(v) The volumetric plastic dissipation leading to reduced heating is primarily due to plastic dilation (expansions). Despite including dilatant plasticity for the mantle, our results did not indicate that volumetric plastic dissipation had significant effects in the hotter parts of the mantle as we did not include assumptions of volumetric plastic strains during ductile creep. Our new model is therefore able to switch between viscoelastic behaviour (controlled by the temperature) and viscoplastic behaviour (controlled by the Drucker–Prager yield criterion which occurs in the brittle regimes).

(vi) We also successfully applied our new model to lithospheric-scale dynamics where we found that whether shear heating and shear with dilatant heating are accounted for, the onset of buckling develops at the same time. However, where we observed a higher volumetric plastic strain, we observed a temperature reduction of
28–33 K, which is about 4 per cent reduction of shear heating when dilatant terms are included.

While the actual temperature change is not so large, we have shown that the impact of accounting for volumetric plastic stresses and strains on several aspects of crustal-scale and lithospheric-scale geodynamic models is not negligible, such as the depths to isotherms, brittle-ductile transition; deformation history, stress evolution and heat flow. Our new model also allowed the emergence of several complexities which were not imposed *a priori*, for example, softening behaviour. The fact that the brittle-ductile transition is a function of the amount of plastic dilatancy also opens new vistas to other questions, for example, the influence on mantle and crustal solidus, thermal steady-state, that is, where the amount of heat produced is evacuated by conduction thereby keeping the temperature in equilibrium; and investigating areas with high heat flows that can be targeted for geothermal exploration.

## ACKNOWLEDGMENTS

This work was supported by Géoscience Environment Toulouse/Observatoire Midi-Pyrénées, Institut Physique du Globe de Paris (IPGP), INTERREG V Caraïbes program, the European Regional Development Fund (FEDER through the PREST project and the AXA Research Fund. EM and ST acknowledge the funding support. HSB appreciates the European Research Council grant PERSISMO (grant number 865411) for partially supporting this work. Dr Nadaya Cubas and Professor Laetitia Le Pourhiet are greatly appreciated for their support and insights; as well as Professor Thibault Duretz for offering useful suggestions which enhanced the quality of the paper. Special appreciation to Dr Antoine Jacquey who thoroughly reviewed an earlier version of this manuscript, and one anonymous reviewer who motivated us to benchmark our new model. Finally, we thank Dr Dàniel Kiss, the two anonymous reviewers of the revised manuscript, the assistant editor, Dr Louise Alexander and the editor, Dr Tobias Keller, who copiously reviewed this work. This is IPGP contribution number 4292.

## DATA AVAILABILITY

All codes developed in the framework of the work and input files to reproduce the experiments reported here are stored in https://github.com/ekemomoh007/VEVP-Code-and-Input-Files.

## REFERENCES

* Abaqus, 2019. Abaqus user subroutines guide, *Dassault Systèmes*, Providence, RI.

* Abaqus, 2010. Abaqus theory manual, *Dassault Systèmes*, Providence, RI.

* Abu Al-Rub, R. & Tehrani, A., 2005. Thermo-mechanical viscoelastic, viscoplastic, and viscodamage model for polymers and polymer composites, *52nd AIAA/ASME/ASCE/AHS/ASC Structures, Structural Dynamics and Materials Conference 19th AIAA/ASME/AHS Adaptive Structures Conference*, 13t (2173).

* Alejano, L.R. & Alonso, E., 2005. Considerations of the dilatancy angle in rocks and rock masses, *Int. J. Rock Mech. Mining Sci.*, **42**(4), 481–507.

* Alejano, L.R. & Bobet, A., 2012. Drucker-Prager criterion, in *The ISRM Suggested Methods for Rock Characterization, Testing and Monitoring: 2007–2014*, pp. 277–296, ed. Ulusay, R., Springer.

* Arthur, J.R.F., Dunstan, T., Al-Ani, Q.A.J. & Assadi, A., 1977. Plastic deformation and failure of granular media, *Géotechnique*, **27**(1), 53–74.

* Babeyko, A.Y. & Sobolev, S.V., 2008. High-resolution numerical modeling of stress distribution in visco-elasto-plastic subducting slabs, *Lithos*, **103**, 205–216.

* Baes, M., Govers, R. & Wortel, R., 2011. Subduction initiation along the inherited weakness zone at the edge of a slab: insights from numerical model, *Geophys. J. Int.*, **184**(3), 991–1008.

* Ben-Israel, A., 1966. A Newton-Raphson method for the solution of systems of equations, *J. Math. Anal. Appl.*, **15**(2), 243–252.

* Bessat, A., Duretz, T., Hetényi, G., Pilet, S. & Schmalholz, S.M., 2020. Stress and deformation mechanisms at a subduction zone: insights from 2-D thermomechanical numerical modelling, *Geophys. J. Int.*, **221**(3), 1605–1625.

* Besson, J., Cailletaud, G., Chaboche, JL., Forest, S., Blètry, M. & Gladwell, G.M.L., 2010. Strain localization phenomena, *Nonlin. Mech. Mater.*, **167**, 371–406.

* Bower, A.F., 2009. *Applied Mechanics of Solids*, CRC Press.

* Boyle, J.T. & Spence, J., 1983. *Stress Analysis for Creep*, Elsevier.

* Branlund, J.M., Regenauer-Lieb, K. & Yuen, D.A., 2000. Fast ductile failure of passive margins from sediment loading, *Geophys. Res. Lett.*, **27**(13), 1989–1992.

* Branlund, J.M., Regenauer-Lieb, K. & Yuen, D.A., 2001. Weak zone formation for initiating subduction from thermo-mechanical feedback of low-temperature plasticity, *Earth planet. Sci. Lett.*, **190**, 237–250.

* Bridgman, P.W., 1964. *Volume changes in the Plastic Stages of Simple Compression, Collected Experimental Volume VII*, pp. 3987–3998, Harvard University Press, Cambridge: MA and London, England.

* Burov, E. & Cloetingh, S., 2010. Plume-like upper mantle instabilities drive subduction initiation, *Geophys. Res. Lett.*, **37**(3), doi: 10.1029/2009GL041535.

* Burov, E., 2011. Rheology and strength of the lithosphere, *Mar. Petrol. Geol.*, **28**(8), 1402–1443.

* Burov, E., Francois, T., Yamato, P. & Wolf, S., 2014. Mechanisms of continental subduction and exhumation of HP and UHP rocks, *Gondw. Res.*, **25–2**, 464–493.

* Capitanio, F., Morra, G. & Goes, S., 2007. Dynamic models of downgoing plate-buoyancy driven subduction: Subduction motions and energy dissipation, *Earth planet. Sci. Lett.*, **262**(1–2), 284–297.

* Chapra, S.C. & Canale, R.P., 2010. *Numerical Methods for Engineers*, McGraw Hill Higher Education.

* Cloetingh, S. A. P. L., Wortel, M. J. R. & Vlaar, N. J., 1982. Evolution of passive continental margins and initiation of subduction zones, *Nature*, **297**(5862), 139–142.

* Cloetingh, S., Wortel, R. & Vlaar, N.J., 1989. On the initiation of subduction zones, in *Subduction Zones Part II*, pp. 7–25, eds, Ruff, L. J. & Kanamori, H., Birkhäuser.

* Cook, N. G. W., 1970. An experiment proving that dilatancy is a pervasive volumetric property of brittle rock loaded to failure, *Rock Mech.*, **2**(4), 181–188.

* Crameri, F. & Kaus, B.J., 2010. Parameters that control lithospheric-scale thermal localization on terrestrial planets, *Geophys. Res. Lett.*, **37**(9), doi: 10.1029/2010GL042921.

* De Borst, R., Crisfield, M. A., Remmers, J. J. & Verhoosel, C. V., 2012. *Nonlinear Finite Element Analysis of Solids and Structures*, John Wiley & Sons.

* de Souza Neto, E.A., 2004. A simple robust numerical integration algorithm for a power-law visco-plastic model under both high and low ratesensitivity, *Commun. Numer. Meth. En.*, **20**(1), 1–17.

* de Souza Neto, E.A., Perić, D. & Owen, D.R.J., 2011. *Computational Methods for Plasticity: Theory and Applications*, Wiley.

* Desai, C.S. & Zhang, D., 1987. Viscoplastic model for geologic materials with generalized flow rule, *Int. J. Numer. Anal. Meth. Geomech.*, **11**, 603–620.

* Drucker, D.C. & Prager, W., 1952. Soil mechanics and plastic analysis or limit design, *Q. Appl. Math.*, **10**, 157–165.

* Duretz, T., Schmalholz, S.M., Podladchikov, Y.Y. & Yuen, D.A., 2014. Physics-controlled thickness of shear zones caused by viscous heating: implications for crustal shear localization, *Geophys. Res. Lett.*, **41**, 4904–4911.
Duretz, T., Schmalholz, S.M. & Podladchikov, Y.Y., 2015. Shear heating-induced strain localization across the scales, *Phil. Mag.*, **95**(28–30), 3192–3207.

Duretz, T., Petri, B., Mohn, G., Schmalholz, S. M., Schenker, F. L. & Müntener, O., 2016. The importance of structural softening for the evolution and architecture of passive margins, *Sci. Rep.*, **6**(1), 38704, doi:10.1038/srep38704.

Duretz, T., Souche, A., de Borst, R. & Le Pourhiet, L., 2018. The benefits of using a consistent tangent operator for viscoelastoplastic computations in geodynamics, *Geochem. Geophys. Geosyst.*, **19**(12), 4904–4924.

Duretz, T., de Borst, R., Yamato, P. & Le Pourhiet, L., 2020. Toward robust and predictive geodynamic modeling: the way forward in frictional plasticity, *Geophys. Res. Lett.*, **47**, e2019GL086027, doi:10.1029/2019GL086027.

Duretz, T., de Borst, R. & Yamato, P., 2021. Modeling lithospheric deformation using a compressible visco-elasto-viscoplastic rheology and the effective viscosity approach, *Geochem. Geophys. Geosyst.*, **22**, e2021GC009675, doi:10.1029/2021GC009675.

Dyksterhuis, S., Albert, R. A. & Müller, R.D., 2005. Finite-element modeling of contemporary and palaeo-intraplate stress using ABAQUS™, *Comput. Geosci.*, **31**(3), 297–307.

England, P.R. & Katz, R., 2010. Melting above the anhydrous solidus controls the location of volcanic arcs, *Nature*, **467**(7316), 700–703.

Ezenwaka, K. C., 2023. The role of subducting heterogeneities and fluids on the Northern Lesser Antilles margin deformation and seismogenesis, *Doctoral dissertation*, Université Côte d'Azur.

Garel, F., Goes, S., Davies, D. R., Davies, J. H., Kramer, S. C. & Wilson, C. R., 2014. Interaction of subducted slabs with the mantle transition-zone: A regime diagram from 2-D thermo-mechanical models with a mobile trench and an overriding plate., *Geochem. Geophys. Geosyst.*, **15**(5), 1739–1765.

Gerbault, M., Poliakov, A. N. & Daignieres, M., 1998. Prediction of faulting from the theories of elasticity and plasticity: what are the limits?, *J. Struct. Geol.*, **20**(2–3), 301–320.

Gerbault, M., 2000. At what stress level is the central Indian Ocean lithosphere buckling?, *Earth Planet. Sci. Lett.*, **178**(1–2), 165–181.

Gerbault, M., Davey, F. & Henrys, S., 2002. Three-dimensional lateral crustal thickening in continental oblique collision: an example from the Southern Alps, New Zealand, *Geochem. J. Int.*, **150**(3), 770–779.

Gerya, T.V. & Yuen, D.A., 2003. Characteristics-based marker-in-cell method with conservative finite-differences schemes for modeling geological flows with strongly variable transport properties, *Phys. Earth planet. Inter.*, **140**(4), 293–318.

Gerya, T.V., Yuen, D.A. & Maresch, W.V., 2004. Thermomechanical modeling of slab detachment, *Earth planet. Sci. Lett.*, **226**(1–2), 101–116.

Gerya, T., 2019. *Introduction to Numerical Geodynamic Modelling*, Cambridge Univ. Press, Cambridge, UK.

Goetze, C. & Evans, B., 1979. Stress and temperature in the bending lithosphere constrained by experimental rock mechanics, *Geophys. J. Int.*, **59**(3), 463–478.

Gurnis, M., 1992. Rapid continental subsidence following the initiation and evolution of subduction, *Science*, **255**(5051), 1556–1558.

Gurnis, M., Hall, C. & Lavier, L., 2004. Evolving force balance during incipient subduction, *Geochem. Geophys. Geosyst*, **5**(7), doi: 10.1029/2003GC000681.

Hassani, R., Jongmans, D. & Chéry, J., 1997. Study of plate deformation and stress in subduction processes using two-dimensional numerical models, *Geophys. Res.: Solid Earth*, **102**(B8), 17 951–17 965.

Hirth, G. & Kohlstedt, D., 2003. Rheology of the upper mantle and the mantle wedge: a view from the experimentalists, *Geophys. Monogr. Ser*, **138**, 83–106.

Holt, A. F., Becker, T. W. & Buffett, B. A., 2015. Trench migration and overriding plate stress in dynamic subduction models, *Geophys. J. Int.*, **201**(1), 172–192.

Jacquey, A.B. & Cacace, M., 2020. Multiphysics modelling of a brittle-ductile lithosphere: 1. explicit visco-elasto-plastic formulation and its numerical implementation, *J. geophys. Res.: Solid Earth*, **125**(1), e2019JB018474, doi: 10.1029/2019JB018474.

Jacquey, A.B., Rattez, H. & Veveakis, M., 2021. Strain localization regularization and patterns formation in rate-dependent plastic materials with multiphysics coupling, *J. Mech. Phys. Solids*, **152**, 104422, doi:10.1016/j.jmps.2021.104422.

Jaquet, Y., Duretz, T. & Schmalholz, S.M., 2016. Dramatic effect of elasticity on thermal softening and strain localization during lithospheric shortening, *Geophys. J. Int.*, **204**, 780–784.

Jourdon, A., Le Pourhiet, L., Petit, C. & Rolland, Y., 2018. Impact of range-parallel sediment transport on 2D thermo-mechanical models of mountain belts: Application to the Kyrgyz Tien Shan, *Terra Nova*, **30**(4), 279–288.

Karrech, A., Regenauer-Lieb, K. & Poulet, T., 2011. Continuum damage mechanics for the lithosphere, *J. geophys. Res.: Solid*, **116**(B4), doi:10.1029/2010JB007501.

Kaus, B.J.P., 2010. Factors that control the angle of shear bands in dynamic numerical models of brittle deformation, *Tectonophysics*, **484**, 36–47.

Kaus, B.J.P. & Podladchikov, Y.Y., 2006. Initiation of localized shear zones in viscoelastoplastic rocks, *J. geophys. Res.*, **111**, B04412, doi:10.1029/2005JB003652.

Kohlstedt, D.L., Evans, B. & Mackwell, S.J., 1995. Strength of the lithosphere: constraints imposed by laboratory experiments, *J. geophys. Res.*, **100**(B9), 17 587–17 601.

Kiss, D., Candidoti, L.G., Duretz, T. & Schmalholz, S.M., 2020. Thermal softening induced subduction initiation at a passive margin, *Geophys. J. Int.*, **220**(3), 2068–2073.

Lallemand, S. & Arcay, D., 2021. Subduction initiation from the earliest stages to self-sustained subduction: insights from the analysis of 70 Cenozoic sites, *Earth Sci. Rev.*, **221**, 103779, doi:10.1016/j.earscirev.2021.103779.

Le Pourhiet, L., 2013. Strain localization due to structural softening during pressure sensitive rate independent yielding, *Bull. Soc. Géol. France*, **184**(4–5), 357–371.

Lemaitre, J. & Chaboche, J.L., 1994. *Mechanics of Solid Materials*, Cambridge University Press.

Leng, W. & Gurnis, M., 2011. Dynamics of subduction initiation with different evolutionary pathways, *Geochem. Geophys. Geosyst.*, **12**(2), doi:10.1029/2011GC003877.

Leng, W. & Gurnis, M., 2015. Subduction initiation at relic arcs, *Geophys. Res. Lett.*, **42**(17), 7014–7021.

Leroy, Y. & Ortiz, M., 1989. Finite element analysis of strain localization in frictional materials, *Int. J. Numer. Anal Methods Geomech.*, **13**, 53–74.

Leroy, Y. & Ortiz, M., 1990. Finite element analysis of transient strain localization phenomena in frictional solids, *Int. J. Numer. Anal Methods Geomech.*, **14**, 93–124.

Manga, M. et al., 2012. Heat flow in the Lesser Antilles island arc and adjacent back arc Grenada basin, *Geochem. Geophys. Geosyst.*, **13**(8), 1–19.

Mase, G.E., 1970. *Continuum Mechanics*, McGraw Hill.

McKenzie, D.M., 1977. The initiation of trenches: a finite amplitude instability, in *Talwani, M. & Pitman III, W.C. Island Arcs, Deep Sea Trenches and Back-arc Basins*, 1, American Geophysical Union 57–61.

Moresi, L., Mühlhaus, H.-B., Lemiale, V. & May, D., 2007. Incompressible viscous formulations for deformation and yielding of the lithosphere, *Geol. Soc. Lond. Special Publ.*, **282**(1), 457–472.

Mueller, S. & Phillips, R. J., 1991. On the initiation of subduction, *J. geophys. Res.: Solid Earth*, **96**(B1), 651–665.

Nikolaeva, K., Gerya, T.V. & Marques, F.O., 2010. Subduction initiation at passive margins: numerical modeling, *J. geophys. Res.*, **15**(B3), doi:10.1029/2009JB006549.

Owen, D.R.J. & Hinton, E., 1980. *Finite Elements in Plasticity: Theory and Practice*, Pineridge Press.

Pajang, S., Le Pourhiet, L. & Cubas, N., 2021. The topographic signature of temperature controlled rheological transitions in an accretionary prism, *Solid Earth Discuss.*, **2021**, 1–27.

Pardoen, B., Seyedi, D.M. & Collin, F., 2015. Shear band modelling in cross-anisotropic rocks, *Int. J. Soils Struct.*, **72**, 63–87.
Patocka, V., Cizkova, H. & Tackley, P.J., 2019. Do elasticity and a free surface affect lithospheric stresses caused by upper-mantle convection?, *Geophys. J. Int.*, **216**(3), 1740–1760.

Perzyna, P., 1986. Fundamental problems in viscoplasticity, *Adv. Appl. Mech.*, **9**, 244–368, doi: 10.1016/S0065-2156(08)70009-7.

Poliakov, A.N., Herrmann, H.J., Podladchikov, Y.Y. & Roux, S., 1994. Fractal plastic shear bands, *Fractals*, **2**(4), 567–581.

Poliakov, A.N. & Herrmann, H.J., 1994. Self-organized criticality of plastic shear bands in rocks, *Geophys. Res. Lett.*, **21**(19), 2143–2146.

Poulet, T. & Veveakis, M., 2016. A viscoplastic approach for pore collapse in saturated soft rocks using redback: an open source parallel simulator for rock mechanics with dissipative feedbacks, *Comput. Geosci.*, **74**, 211–221.

Ranalli, G., 1995. *Rheology of the Earth*, Chapman & Business Media.

Ravichandran, G., Rosakis, A.J., Hodowany, J. & Rosakis, P., 2002. On the conversion of plastic work into heat during high-strain-rate deformation, in AIP conference proceedings, *Am. Inst. Phys.*, **620**(1), 557–562.

Regenauer-Lieb, K. & Yuen, D.A., 1998. Rapid conversion of elastic energy into plastic shear heating during incipient necking of the lithosphere, *Geophys. Res. Lett.*, **25**(14), 2737–2340.

Regenauer-Lieb, K., Yuen, D.A. & Branlund, J., 2001. The initiation of subduction: criticality by addition of water, *Science*, **294**(14), 578–234.

Regenauer-Lieb, K. & Yuen, D.A., 2003. Modelling shear zones in geological and planetary sciences: solid- and fluid-thermal-mechanical approaches, *Earth Planet. Sci. Lett.*, **63**, 295–349.

Rittel, D., 1999. On the conversion of plastic work to heat during high strain rate deformation of glassy polymers, *Mech. Mater.*, **31**(2), 131–139.

Rosakis, P., Rosakis, A., Ravichandran, G. & Hodowany, G., 2000. A thermodynamic internal variable model for the partition of plastic work into heat and stored energy in metals, *J. Mech. Phys. Solids*, **48**(3), 581–607.

Ruh, J.B., Le Pourhiet, L., Burov, E. & Gerya, T., 2015. Tectonic slicing of subducting oceanic crust along plate interfaces: numerical modeling, *Geochem. Geophys. Geosyst.*, **16**(10), 3505–3531.

Salomon, C., 2018. Finite element modelling of the geodynamic processes of the Central Andes subduction zone: a reference model, *Geod. Geodyn.*, **9**(3), 246–251.

Schmalholz, S.M., Duretz, T., Schenker, F.L. & Podladchikov, Y.Y., 2014. Kinematics and dynamics of tectonic nappes: 2D numerical modelling and implications for high and ultra-high pressure tectonism in the Western Alps, *Tectonophysics*, **631**, 160–175.

Schmeling, H. et al., 2008. A benchmark comparison of spontaneous subduction models—Towards a free surface, *Phys. Earth planet. Inter.*, **171**(1–4), 198–223.

Simo, J.C. & Hughes, T.J.R., 1998. *Computational Inelasticity*, Springer.

Stern, R. & Gerya, T., 2018. Subduction initiation in nature and models: a review, *Tectonophysics*, **746**, 173–198.

Stern, R.J., 2004. Subduction initiation: spontaneous and induced, *Earth planet. Sci. Lett.*, **226**, 275–292.

Skrzypek, J., 1993. *Plasticity and Creep: Theory, Examples and Problems*, CRC Press.

Thielmann, M. & Kaus, B.J.P., 2012. Shear heating induced lithospheric-scale localization: Does it result in subduction?, *Earth planet. Sci. Lett.*, **359**, 1–13.

Thielmann, M., Rozel, A., Kaus, B.J.P. & Ricard, Y., 2015. Intermediate-depth earthquake generation and shear zone formation caused by grain size reduction and shear heating, *Geology*, **43**, 791–794.

Toth, J. & Gurnis, M., 1998. Dynamics of subduction initiation at preexisting fault zones, *J. geophys. Res.: Solid Earth*, **103**(B8), 18 053–18 067.

Turcotte, D.L. & Schubert, G., 2002. *Geodynamics*, Cambridge Univ. Press.

Vermeer, P., 1990. The orientation of shear bands in biaxial tests, *Géotechnique*, **40**(2), 223–236, doi: 10.1680/geot.1990.40.2.223.

Vermeer, P.A. & de Borst, R., 1984. Non-associated plasticity for soils, concrete and rock, *Heron*, **29**(3), 1–62, Available at: https://typeset.io/pdf/non-associated-plasticity-for-soils-concrete-and-rock-1k1jstx009.pdf.

Vogt, K., Matenco, L. & Cloetingh, S., 2017. Crustal mechanics control the geometry of mountain belts. Insights from numerical modelling, *Earth planet. Sci. Lett.*, **460**, 12–21.

Willis, K., & Hooper, A., 2019. Strain localization by shear heating and the development of lithospheric shear zones, *Tectonophysics*, **764**, 62–76.

Yqma, T.J., 1995. Historical development of the Newton-Raphson method, *SIAM Rev.*, **37**(4), 531–551.

Zhang, L., Zlotnik, S. & Li, C.F., 2021. Anomalous subduction initiation: young under old oceanic lithosphere, *Geochem. Geophys. Geosyst.*, **22**(6), e2020GC009549, doi:10.1029/2020GC009549.

Zhao, X. G. & Cai, M., 2010. A mobilized dilation angle model for rocks, *Int. J. Rock Mech. Mining Sci.*, **47**(3), 368–384.

Zhou, X. & Wada, I., 2021. Differentiating induced versus spontaneous subduction initiation using thermomechanical models and metamorphic soles, *Nat. Commun.*, **12**(1), 1–10.

Zienkiewicz, O.C. & Cormeau, I.C., 1974. Viscoplasticity-plasticity and creep in elastic solids-a unified numerical solution approach, *Int. J. Numer. Methods Eng.*, **8**, 821–845.

Zienkiewicz, O.C. & Taylor, R.L., 2005. *The Finite Element Method for Solid and Structural Mechanics*, Elsevier.

## APPENDIX A: NUMERICAL IMPLEMENTATION OF CONSTITUTIVE LAWS

The formulations to solve the mechanical and thermal diffusion parts of the system described in Section 2 have been coded into a User-defined mechanical material behaviour (UMAT) and user-defined thermal material behaviour UMATHT subroutines, respectively, and implemented in Abaqus® (Abaqus 2010, 2019).

Given a displacement at a given time step which can result from a velocity or force-driven boundary condition, Abaqus® supplies the total strain increment at that time step. Using our proposed constitutive laws, corresponding elastic stresses are computed. The viscous and viscoplastic corrections to the stresses are at the heart of the stress update scheme, and for an adaptive time-stepping scheme that we adopt, the consistent algorithmic tangent modulus is required. We describe the details of our algorithmic implementation here.

### A1 Stress update scheme

At each time step, we first assume that the strain increment, $\Delta\epsilon_{ij}^e$, at that time step is fully elastic. We then compute an elastic trial stress.

$$ \sigma_{ij}^{e, \text{trial}} = C_{ijkl}^{e} (\epsilon_{kl}^{e, \text{trial}} + \Delta\epsilon_{kl}^{e, \text{trial}}). \quad (A1) $$

Given the elastic trial stress, we obtain trial deviatoric stresses and pressures from:

$$ \begin{align} s_{ij}^{e, \text{trial}} &= \sigma_{ij}^{e, \text{trial}} - \frac{\sigma_{kk}^{e, \text{trial}}}{3} \delta_{ij}, \tag{A2} \\ P^{e, \text{trial}} &= \frac{I_{1}^{e, \text{trial}}}{3} = \frac{\sigma_{kk}^{e, \text{trial}}}{3}. \notag \end{align} $$

The trial deviatoric stresses are used to compute the invariant of deviatoric stresses in the elastic state. Given the function in eq. (9), we can write a viscous function in the form:

$$ \Phi^{\nu}(\sigma_{ij}^{e}, T) = (\Phi^{\nu})^{m} f(T) - \frac{\Delta\gamma^{\nu}}{\Delta t} = 0, \quad (A3) $$

where $\Delta t$ is the time increment at a step, $\Phi^{\nu} = \sqrt{J_{\text{II}}^{e}(\sigma_{ij}^{e, \text{trial}})} - G\Delta\gamma^{\nu}$. Eq. (A3) is a nonlinear equation and we obtain the unknown $\Delta\gamma^{\nu}$ by performing Newton-Raphson iterations (Ben-Israel 1966; Yqma 1995) using the derivative of $\Phi^{\nu}$
with respect to $\Delta\gamma^v$:

$$ \frac{d\tilde{\Phi}^v}{d\Delta\gamma^v} = -mG(\Phi^v)^{m-1} f(T)\Delta t - 1 = 0. \quad (A4) $$

After computing $\Delta\gamma^v$, update the deviatoric stresses as follows:

$$ \Delta\epsilon_{ij}^v = \Delta\gamma^v(\sigma_{ij}^{e, \text{trial}}, T) \frac{\partial \Phi_F^v}{\partial s_{ij}^{e, \text{trial}}} = \Delta\gamma^v \frac{s_{ij}^{e, \text{trial}}}{2\sqrt{J_{\text{II}}^e(s_{ij}^{e, \text{trial}})}}, \quad (A5) $$

$$ s_{ij}^v = \left( 1 - \frac{G \Delta \gamma^v}{\sqrt{J_{\text{II}}^e(s_{ij}^{e, \text{trial}})}} \right) s_{ij}^{e, \text{trial}}. \quad (A6) $$

A deviatoric stress invariant ($J_{\text{II}}^v$) is computed after correcting for creep. This new stress invariant is fed to assess the viscoplastic routine through the Drucker-Prager yield function in eq. (11). If the yield criterion is violated, we perform more Newton-Raphson iterations to obtain the viscoplastic multiplier. We therefore write a residual viscoplastic potential function using the functional form of (Zienkiewicz & Cormeau. 1974; Perzyna 1986):

$$ \tilde{\Phi}_{\text{R}}^{\text{vp}}(\sigma_{ij}^{v}, T) = \frac{\Delta t}{\mu} \left( \frac{\Phi_{\text{R}}^{\text{DP}}}{\Phi_{0}} \right)^{m} - \Delta \gamma^{\text{vp}} = 0, \quad (A7) $$

where $\Phi_{\text{R}}^{\text{DP}} = \Phi_{\text{Y}}^{\text{DP}} - c_1 \Delta \gamma^{\text{vp}}$, with $c_1 = (G + \alpha_1 \alpha_3 K)$, $\Phi_0 = c_0$ and $\mu$ is a viscosity-related term. Because eq. (A7) is a nonlinear equation, we find the unknown $\Delta\gamma^{\text{vp}}$ by Newton-Raphson iterations as the creep case. The derivative of the residual viscoplastic potential function required for these iterations is given by:

$$ \frac{d\tilde{\Phi}_{\text{R}}^{\text{vp}}}{d\Delta\gamma^{\text{vp}}} = -\frac{c_1 m \Delta t}{c_0 \mu} \left( \frac{\Phi_{\text{Y}} - c_1 \Delta\gamma^{\text{vp}}}{c_0} \right)^{m-1} - 1. \quad (A8) $$

Eqs. (A7) and (A8) are used to perform Newton-Raphson iterations to obtain the unknown $\Delta\gamma^{\text{vp}}$. A critical aspect of the convergence of the Newton-Raphson iterations for obtaining both the viscous and viscoplastic multipliers is the requirement that the multipliers are strictly positive and that the algorithm converges. To aid the convergence of the Newton-Raphson scheme, we have included a modified bisection method when these conditions are violated (de Souza Neto 2004; Chapra & Canale 2011).

After obtaining the viscoplastic multiplier upon yielding, the stresses, viscoplastic strains, temperature change and state variables can be computed:

$$ \Delta\epsilon_{ij}^{\text{vp}} = \Delta\gamma^{\text{vp}} \left( \frac{s_{ij}^{v}}{2\sqrt{J_{\text{II}}^{v}(s_{ij}^{v})}} + \frac{\alpha_3}{3} \delta_{ij} \right) \quad (A9) $$

$$ s_{ij}^{\text{vp}} = \left( 1 - \frac{G \Delta \gamma^{\text{vp}}}{\sqrt{J_{\text{II}}^{v}(s_{ij}^{v})}} \right) s_{ij}^{v}, \quad (A10) $$

$$ P = P^{\text{e, trial}} - \alpha_3 \Delta \gamma^{\text{vp}} K, \quad (A11) $$

$$ \sigma_{ij} = s_{ij}^{\text{vp}} + P, \quad (A12) $$

$$ \Delta T = \frac{\sigma_{ij} \Delta \epsilon_{ij}^{\text{vp}}}{\rho C_p}. \quad (A13) $$

## A2 Consistent algorithmic tangent modulus (CATM)

Finally, using the Newton scheme to solve the global equations, an important consideration is the convergence of the solver and the selection of an adaptive time step. To this end, we estimate a convergence matrix obtained by linearizing the stress update

schemes for deviatoric stresses and pressure given in eqs (A6), (A10) and (A11). We proceed by taking small increments in the deviatoric stress and pressure with respect to the trial elastic state. We compute this matrix for the creep and viscoplastic states.

### A2.1 Consistent algorithmic tangent modulus for creep deformation

$$ C_{ijkl}^{v} = \frac{d\sigma_{ij}}{d\epsilon_{kl}^{e, \text{trial}}} = \frac{ds_{ij}^{v}}{d\epsilon_{kl}^{e, \text{trial}}} + \delta_{ij} \left( \frac{dP}{d\epsilon_{kl}^{e, \text{trial}}} \right) \quad (A14) $$

The expression for deviatoric stresses can be equivalently written in terms of the deviatoric strains ($e_{ij}^{e, trial}$) taking into account that $s_{ij}^{e, trial} = 2Ge_{ij}^{e, trial}$, and $e_{ij}^{trial} = \epsilon_{ij}^{trial} - \epsilon_{kk}^{trial}/3$:

$$ s_{ij}^{v} = 2G \left( 1 - \frac{\Delta\gamma^{v}}{\sqrt{2}e_{\text{norm}}^{\text{trial}}} \right) e_{ij}^{\text{trial}}. \quad (A15) $$

$e_{\text{norm}}^{\text{trial}} = \sqrt{e_{ij}^{\text{trial}} e_{ij}^{\text{trial}}}$ represents the Euclidean norm of the deviatoric strain tensor. The elements to assemble the consistent Jacobian matrix during creep deformation are given by derivatives of the stress states after creep:

$$ \frac{ds_{ij}^{v}}{d\epsilon_{kl}^{e, \text{trial}}} = 2G \left( 1 - \frac{\Delta\gamma^{v}}{\sqrt{2}e_{\text{norm}}^{\text{trial}}} \right) I_{ijkl}^{\text{dev}} + \sqrt{2}G \left( \frac{\Delta\gamma^{v}}{e_{\text{norm}}^{\text{trial}}} - b_1 \right) \frac{e_{ij}^{\text{trial}} e_{kl}^{\text{trial}}}{e_{\text{norm}}^{\text{trial}} e_{\text{norm}}^{\text{trial}}}, \quad (A16) $$

$$ \delta_{ij} \frac{dP^{e, trial}}{d\epsilon_{kl}^{e, trial}} = K \delta_{ij} \delta_{kl}, \quad (A17) $$

with

$$ b_1 = \frac{\sqrt{2G (\Delta tf(T))^{1/m}}}{\frac{(\Delta \gamma^v)^{(1-m)/m}}{m} + G(\Delta tf(T))^{1/m}}. $$

The expression for the consistent Jacobian, when creep deformation is considered without viscoplasticity, is given by a combination of Equations (A16) and (A17):

$$ C_{ijkl}^{v} = 2G \left( 1 - \frac{\Delta \gamma^{v}}{\sqrt{2} e_{\text{norm}}^{\text{trial}}} \right) I_{ijkl}^{\text{dev}} + \sqrt{2} G \left( \frac{\Delta \gamma^{v}}{e_{\text{norm}}^{\text{trial}} - b_1} \right) \frac{e_{ij}^{\text{trial}} e_{kl}^{\text{trial}}}{e_{\text{norm}}^{\text{trial}} e_{\text{norm}}^{\text{trial}}} + K \delta_{ij} \delta_{kl}, \quad (A18) $$

where $I_{ijkl}^{\text{dev}} = (\delta_{ik}\delta_{jl} + \delta_{il}\delta_{jk})/2 - \delta_{ij}\delta_{kl}/3$ is the fourth order deviatoric projection tensor, which extracts the deviatoric components of a stress or strain tensor.

### A2.2 Consistent algorithmic tangent modulus for viscoplastic deformation

The linearization scheme follows the same procedure used for the creep state. For the deviatoric stresses,

$$ \frac{ds_{ij}^{\text{vp}}}{d\epsilon_{kl}^{e, \text{trial}}} = \left( 1 - \frac{G \Delta \gamma^{\text{vp}}}{\sqrt{J_{\text{II}}^{v}(s_{ij}^{v})}} \right) \frac{ds_{ij}^{v}}{d\epsilon_{kl}^{e, \text{trial}}}. \quad (A19) $$
$$ \frac{ds_{ij}^{vp}}{ds_{kl}^{e \text{trial}}} = \left(1 - \frac{G \Delta \gamma^{vp}}{\sqrt{J_{H}^{v}(s_{ij}^{v})}}\right) \frac{ds_{ij}^{v}}{ds_{kl}^{e \text{trial}}} - \frac{\alpha_{1} b_{2} K G}{\sqrt{J_{H}^{v}(s_{ij}^{v})}} \delta_{ij} \frac{s_{kl}^{v}}{\sqrt{J_{H}^{v}(s_{kl}^{v})}} - \frac{G}{\sqrt{J_{H}^{v}(s_{ij}^{v})}} \left(\sqrt{2} G b_{2} \left(1 - \frac{b_{1}}{\sqrt{2}}\right) \frac{e_{ij}^{\text{trial}}}{e_{\text{norm}}} + \frac{\sqrt{2} G \Delta \gamma^{vp}}{\sqrt{J_{H}^{v}(s_{ij}^{v})}} \left(\frac{b_{1}}{\sqrt{2}} - 1\right) \frac{e_{ij}^{\text{trial}}}{e_{\text{norm}}}\right) \frac{s_{kl}^{v}}{\sqrt{J_{H}^{v}(s_{kl}^{v})}} \\ \frac{\delta_{ij} \frac{d P_{n+1}}{d s_{kl}^{e \text{trial}}}}{ds_{kl}^{e \text{trial}}} = K (1 - \alpha_{1} \alpha_{3} b_{2} K) \delta_{ij} \delta_{kl} + \alpha_{3} b_{2} K (b_{1} - \sqrt{2} G) \frac{e_{ij}^{\text{trial}}}{e_{\text{norm}}} \delta_{kl}. \quad (A20) $$

where

$$ b_2 = \frac{\frac{1}{c_0} \left( \frac{\Delta t}{\mu} \right)^{1/m}}{\frac{(\Delta \gamma^{vp})^{(1-m)/m}}{m} + \frac{c_1}{c_0} \left( \frac{\Delta t}{\mu} \right)^{1/m}} $$

The procedure of obtaining the tangent moduli is a linearization of the stress update scheme, such that the adaptive time stepping is a function of the deformation state of the material.