# A Comparison of Plasticity Regularization Approaches for Geodynamic Modeling

Thibault Duretz, Ludovic Räss, René de Borst, Tim Hageman

► To cite this version:

Thibault Duretz, Ludovic Räss, René de Borst, Tim Hageman. A Comparison of Plasticity Regularization Approaches for Geodynamic Modeling. Geochemistry, Geophysics, Geosystems, 2023, 24 (7), pp.e2022GC010675. 10.1029/2022GC010675. insu-04181518

HAL Id: insu-04181518

https://insu.hal.science/insu-04181518

Submitted on 16 Aug 2023

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
# Geochemistry, Geophysics, Geosystems®

## RESEARCH ARTICLE

10.1029/2022GC010675

**Special Section:**

Frontiers in lithospheric dynamics: bridging scales through observations, experiments, and computations

**Key Points:**

* Regularization approaches for plastic strain localization are tested using a single code based on pseudo-transient method

* All considered schemes provide convergent result upon mesh refinement and satisfaction of equilibrium

* The use of Cosserat continuum is most robust regularization approach and is also is the most demanding

**Correspondence to:**

T. Duretz,
duretz@uni-frankfurt.de

**Citation:**

Duretz, T., Räss, L., de Borst, R., & Hageman, T. (2023). A comparison of plasticity regularization approaches for geodynamic modeling. *Geochemistry, Geophysics, Geosystems, 24*, e2022GC010675. https://doi.org/10.1029/2022GC010675

Received 30 AUG 2022
Accepted 20 MAY 2023

© 2023 The Authors. Geochemistry, Geophysics, Geosystems published by Wiley Periodicals LLC on behalf of American Geophysical Union. This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

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

## A Comparison of Plasticity Regularization Approaches for Geodynamic Modeling

Thibault Duretz¹,² ID, Ludovic Räss³,⁴ ID, René de Borst⁵ ID, and Tim Hageman⁶ ID

¹Institut für Geowissenschaften, Goethe-Universität Frankfurt, Frankfurt, Germany, ²CNRS, Géosciences Rennes UMR 6118, University Rennes, Rennes, France, ³Laboratory of Hydraulics, Hydrology and Glaciology (VAW), ETH Zurich, Zurich, Switzerland, ⁴Swiss Federal Institute for Forest, Snow and Landscape Research (WSL), Birmensdorf, Switzerland, ⁵Department of Civil and Structural Engineering, University of Sheffield, Sheffield, UK, ⁶Department of Civil and Environmental Engineering, Imperial College London, London, UK

**Abstract** The emergence, geometry and activation of faults are intrinsically linked to frictional rheology. The latter is thus a central element in geodynamic simulations which aim at modeling the generation and evolution of fault zones and plate boundaries. However, resolving frictional strain localization in geodynamic models is problematic. In simulations, equilibrium cannot always be attained and results can depend on mesh resolution. Spatial and temporal regularization techniques have been developed to alleviate these issues. Herein, we investigate three popular regularization techniques, namely viscoplasticity, gradient plasticity and the use of a Cosserat continuum. These techniques have been implemented in a single framework based on an accelerated pseudo-transient solution strategy. The latter allows to explore the effects of regularization on shear banding using the same code and model configuration. We have used model configurations that involve three levels of complexity: from the emergence of a single isolated shear band to the visco-elasto-plastic stress buildup of a crust. All considered approaches allow to resolve shear banding, provide convergence upon mesh refinement and satisfaction of equilibrium. Viscoplastic regularization is straightforward to implement in geodynamic codes. Nevertheless, more stable shear banding patterns and strength estimates are achieved with computationally more expensive gradient and Cosserat-type regularizations. We discuss the relative benefits of these techniques and their combinations for geodynamic modeling. Emphasis is put on the potential of Cosserat-type media for geodynamic applications.

### 1. Introduction

Strain localization is ubiquitous and controls the development of shear zones and the establishment of plate boundaries. The outermost portion of the lithosphere, the crust, is relatively cold. In these regions, deformation is governed by elasto-plasticity, which allows tectonic plates to behave both as rigid solid and, locally, as a deforma-ble solid. The rocks that make up the crust exhibit a frictional behavior, which causes their strength to depend on pressure (Byerlee, 1978). Moreover, the dilatancy angles of rocks are generally smaller than friction angles (Zhao & Cai, 2010), therefore frictional plasticity models should be non-associated. Geodynamic models are typically employed to simulate the deformations of the Earth's upper shell, such as the opening of rifts zone (e.g., Naliboff et al., 2017), the growth of orogenic wedges (e.g., Buiter et al., 2016), or the formation of transform boundaries (e.g., Gerya, 2013) heavily rely on the simulation of frictional deformation processes. Nonetheless, the inclusion of frictional plasticity in numerical models poses problems as standard model implementations may not achieve force equilibrium and may exhibit mesh dependence. The absence of characteristic spatial or temporal scales in standard frictional constitutive laws is the primary cause for these problems. To overcome these issues, models should be augmented by including additional physics (e.g., fluid pressure diffusion, viscoplastic dissipation) (e.g., Brantut et al., 2017), which from a mathematical point of view can be considered as regularization tech-niques (e.g., de Borst et al., 1993).

Numerous regularization strategies have been proposed, which have the effect of introducing a length or a time scale in the model. Viscoplastic regularization relies on the inclusion of a time scale (de Borst & Duretz, 2020; W. M. Wang et al., 1997) which implicitly rather than directly introduces a length scale like in (e.g., H. Wang, 2019). The yield function does not depend on the gradients of plastic strain: it is defined locally. As a consequence, return mapping can be achieved without solving an additional partial differential equation. Thus, the total number of degrees of freedom is not increased. The implementation of viscoplastic regularization in existing codes is
straightforward and the approach has been used in the context of geodynamic modeling (Duretz et al., 2019, 2020; Jacquey & Cacace, 2020). However, during dynamic deformation viscoplasticity alone is insufficient to remove mesh sensitivity (Stathas & Stefanou, 2021). Hence, there is a need to investigate the relative benefits of other, in particular spatial regularization techniques for modeling shear banding.

Non-local plasticity is one of the earliest regularization approaches (Bažant & Lin, 1988). It involves the spatial averaging of the plastic strain, which is used for in hardening/softening laws. A characteristic length scale, proportional to the area over which strain is averaged is introduced. With this scheme, the sparsity of the linearized discrete operators is reduced as the length scale is increased. Since this scheme is based on strain averaging, it is a candidate for regularizing mesh sensitivity which arises from material strain softening. Gradient plasticity involves the spatial gradients of the plastic strain, typically second-order gradients (i.e., Laplacian), in the yield function (de Borst & Mühlhaus, 1992). This implies the existence of a multiplier ($K_g$) in unit Pa·m² and introduces a length scale in the constitutive model. Due to this modification, the return-mapping procedure involves the solution of a partial differential equation and the plastic multiplier rate becomes an additional degree of freedom (DOF).

Spatial regularization can also achieved by employing a Cosserat, or more generally a micro-polar continuum. This approach relies on the inclusion of micro-rotations and can capture the micro-structure of granular or blocky materials. The constitutive relation which relates micro-curvatures and the conjugate couple-stresses introduces a length scale (Mühlhaus & Vardoulakis, 1987; Sabet & de Borst, 2019; Stefanou et al., 2017). Different from viscoplasticity and gradient plasticity, the use of a Cosserat continuum for regularization does not require modification of the yield function, but the second and third stress invariants have to be changed. In contrast, tensor invariants are augmented by the including contributions from micro-curvatures and couple-stresses. Another important difference is the fact that, when using a Cosserat continuum, length scale effects do not only affect the plastic strains, but also their elastic and viscous counterparts. However, the use of a Cosserat continuum is only effective for shear localization, since the rotations are mobilized only in mode-II.

The use of gradient and Cosserat plasticity has so far mostly been limited to engineering and geomechanical modeling (de Borst et al., 1993; Mühlhaus & Vardoulakis, 1987; Sabet and de Borst, 2019; Stefanou et al., 2017). The benefits and drawbacks of each approach have been discussed in independent studies, which involve different model configurations and different simulation tools. While being insightful, comparison of regularization techniques using similar tools and similar model configurations are rare (de Borst et al., 1993).

Here, we investigate the effects of the three main regularization techniques for non-associated frictional plastic deformation (viscoplasticity, second-order gradient, and Cosserat plasticity) in a similar context. We use the finite difference method (FDM) and an accelerated pseudo-transient (PT) solution strategy (Räss et al., 2022), which provides a simple and unified framework to study coupled non-linear systems (Duretz, Räss, et al., 2018; Räss et al., 2019). We also investigate the combination of spatial and temporal regularization approaches. This has been applied in damage mechanics, for example, through the phase-field approach (e.g., Miehe et al., 2010), and has also been proposed for the study of rate-and-state frictional slip (e.g., Pranger et al., 2022) and compaction banding (e.g., Leuthold et al., 2021).

## 2. Model Description

### 2.1. Governing Equations

We consider steady-state deformations of a compressible medium in two-dimensions (2D plane strain). The balance of linear momentum takes the form of:

$$ \begin{gathered} \frac{\partial \sigma_{xx}}{\partial x} + \frac{\partial \sigma_{yx}}{\partial y} + \rho g_x = 0 \\ \frac{\partial \sigma_{xy}}{\partial x} + \frac{\partial \sigma_{yy}}{\partial y} + \rho g_y = 0 \end{gathered} \qquad (1) $$

where body force components are products of the density ($\rho$) and the gravity acceleration vector components ($g_i$). The total stress is expressed as $\sigma_{ij} = -p\delta_{ij} + \tau_{ij} + \epsilon_{ijk}R_k$ which contains contributions from the symmetric part of the deviatoric stress tensor ($\tau_{ij}$), the pressure ($p$) and the force-conjugate of the micro-rotation rate $\dot{\omega}_k$. The latter contributions is antisymmetric contribution and is non-zero only in the case of the Cosserat model.
The symbol $c_{ijk}$ represents the Levi-Civita tensor. In the Cosserat model, since the stress tensor is not symmetric, the balance of angular momentum also needs to be considered:

$$ \frac{\partial m_{xz}}{\partial x} + \frac{\partial m_{yz}}{\partial y} - (\sigma_{yx} - \sigma_{xy}) = 0 \qquad (2) $$

where $m_{ij}$ are components of the couple-stress tensor.

The conservation of mass is formulated as

$$ \left( \frac{\partial v_x}{\partial x} + \frac{\partial v_y}{\partial y} \right) = \dot{\epsilon}_{kk}^e + \dot{\epsilon}_{kk}^p \qquad (3) $$

where $v_i$ are velocity vector components. The terms $\dot{\epsilon}_{kk}^e$ and $\dot{\epsilon}_{kk}^p$ correspond to elastic and plastic divergence rates, respectively. In the following, we will solve Equations 1–3 to obtain the velocity vector components, the pressure and the micro-rotation-rate vector $\dot{\omega}_i$, given initial and boundary conditions, as well as the constitutive relationships described below. In the following computations, we will consider fixed normal velocity components to model boundaries, zero shear rate tangential to boundaries and zero micro-rotation rate components at the boundaries.

## 2.2. Constitutive Relationships

In order to capture long-term creep of rocks, strain localization due to frictional plastic deformation and rock elasticity, our rheological model involves contributions from viscosity ($\nu$), elasticity ($e$) and plasticity ($p$). The deviatoric rheology is based on an additive decomposition of the deviatoric strain rate tensor ($\dot{\epsilon}_{ij}$):

$$ \dot{\epsilon}_{ij} = \dot{\epsilon}_{ij}^v + \dot{\epsilon}_{ij}^e + \dot{\epsilon}_{ij}^p = \frac{\tau_{ij}}{2\eta} + \frac{\dot{\tau}_{ij}}{2G} + \dot{\lambda} \frac{\partial Q}{\partial \tau_{ij}}, \qquad (4) $$

where $\eta$ is the creep viscosity, $G$ is the shear modulus, $\dot{\lambda}$ is the rate of the plastic multiplier and $Q$ is the plastic flow potential. For the volumetric rheology, we consider elastic and plastic deformations:

$$ \dot{\epsilon}_{kk} = \dot{\epsilon}_{kk}^e + \dot{\epsilon}_{kk}^p = -\frac{\dot{p}}{K} - \dot{\lambda} \frac{\partial Q}{\partial p}, \qquad (5) $$

where $K$ is the bulk modulus. For Cosserat continuum, the following decompositions are applied as well:

$$ \dot{\kappa}_{ij} = \dot{\kappa}_{ij}^v + \dot{\kappa}_{ij}^e + \dot{\kappa}_{ij}^p = \frac{m_{ij}}{2\eta l_c^2} + \frac{\dot{m}_{ij}}{2G l_c^2} + \dot{\lambda} \frac{\partial Q}{\partial m_{ij}} \qquad (6) $$

$$ \dot{W}_k = \dot{W}_k^v + \dot{W}_k^e + \dot{W}_k^p = \frac{R_k}{2\eta_c} + \frac{\dot{R}_k}{2G_c} + \dot{\lambda} \frac{\partial Q}{\partial R_k}, \qquad (7) $$

where $\dot{\kappa}_{ij}$ are the curvature rates, $l_c$ is the Cosserat length scale, and $\eta_c$ and $G_c$ are the Cosserat viscosity and the shear modulus, respectively. $\dot{W}_k$ is the rotation rate which contains contributions from the spin ($\dot{\omega}_k$) and micro-rotation rate ($\dot{\omega}_k$). The kinematic and semi-discrete constitutive relationships are explicitly given in Appendices A and B.

We use a Drucker-Prager model, and the yield function and the plastic potential are expressed as, respectively:

$$ F = \tau_{II} - p \sin \psi - C \cos \phi, \qquad (8) $$

$$ Q = \tau_{II} - p \sin \psi, $$

where $\phi$ is the friction angle, $C$ is the cohesion, $\psi$ is the dilatancy angle and $\tau_{II}$ is the square root of the second invariant of the deviatoric stress tensor. The progressive strain-induced decohesion is modeled by applying linear strain softening:

$$ \dot{C} = \sqrt{\frac{2}{3}} h \dot{\lambda}, \qquad (9) $$
### 3. Pseudo-Transient (PT) Solving Strategy

The above equations were discretized with the FDM using a staggered-grid spatial discretization and a backward-Euler temporal discretization. The resulting set of non-linear equations was solved using the accelerated PT method (Räss et al., 2022) which relies on an explicit second-order integration of the non-linear equations in pseudo-time. This approach is particularly well-suited for the solution of coupled non-linear equations (Duretz, Räss, et al., 2018; Räss et al., 2017, 2022) and its flexibility allows to easily incorporate additional equations, combine new elements and explore the effects of coupled physical processes (e.g., Schmalholz et al., 2020). In an accelerated PT solve, damping is applied to each equation that include a Laplace-type operator (Räss et al., 2022). In the present case, damping was thus applied to the linear momentum balance, the angular momentum balance (for Cosserat continuum) and the return-mapping equation (for gradient regularization). Besides damping, a continuation was applied to the plastic multiplier. This allows to progressively relax the non-linearity due to plasticity throughout the iterative process in a similarly fashion to treatment of power-law viscosity in, for example, Duretz, Räss, et al. (2018). The values of all involved numerical parameters are given in the appended code.

We implement the discretized governing equations using the Julia language (Bezanson et al., 2017). We use the ParallelStencil.jl (Omlin & Räss, 2021b) and ImplicitGlobalGrid.jl (Omlin & Räss, 2021a) Julia packages which permit to write backend-agnostic high-level code for high-performance distributed stencil computations on xPUs, that is, both central processing units and graphical processing units (GPUs). This approach conveniently addresses the *two-language problem* allowing for the development of a single code that can be used both for prototyping and production purpose. The accelerated PT method provides a fully local and iterative algorithm optimally leverages the processing capabilities of many-core hardware such as GPUs.

### 4. Model Configuration

Next, we investigate shear banding using three different model configurations. The first case focuses on the development of a single shear band. A 14.1 × 10.4 km² model domain containing a weak seed in the southwest corner (radius *R* and shear modulus *G*/4) is subjected to background pure shear rate ($\dot{\epsilon}^{BG}$). The shear stress and the micro-rotation rate are zero on the domain boundary. The gradient of the plastic strain at the boundary (i.e., variable in space and time) is set to zero. We do not aim at resolving the complete elastic loading, and for this reason the models are pre-stressed. The initial pressure is set to the confining pressure ($P_{conf} = 250$ MPa), the horizontal deviatoric stress is set to $-P_{conf}/2$, and the horizontal deviatoric stress is set to $P_{conf}/2$. For models with a Cosserat continuum, the initial couple-stress components and bending stress are set to zero.

The second case involves a domain size of 14.0 × 10 km². The boundary conditions are similar to the previous case. The initial condition accounts for a random cohesion field which gives rise to a non-trivial shear banding pattern. All models were run with the same initial random noise, which was interpolated to meshes of different resolution. The initial cohesion is set to 100 MPa. In order to introduce a spatial variation, a read noise is applied to the cohesion field. The latter is isotropic and is characterized by a zero mean, a constant variance and no preferred period. Diffusion is applied to the initial cohesion for a total time of 7.1 Kyr and with a diffusivity of 10⁻⁶ m² · s⁻¹ resulting in minimum and maximum cohesion values of 95.6 and 103.8 MPa, respectively. In crustal-scale models, the domain has dimensions of 50 × 30 km². The vertical component acceleration of gravitational acceleration ($g_y$) and density are respectively set to -9.81 m/s² and 2,700 kg/m³. The initial pressure is set to lithostatic pressure ($P = -\int \rho g_y dy$) and the top surface is a free surface. All initial stresses are set to zero. The temperature (*T*) increases with depth with a gradient of 15°C/km and the surface temperature is 20°C. Viscous creep is activated and the viscosity follows a power-law relation, $\eta^v = 2^{-n} A \exp(-Q_a/RT) \tau_{II}^n$, using creep parameters of Westerly granite (Hansen and Carter, 1983): $A = 3.1623 \times 10^{-26}$ Pa⁻⁵.s⁻¹, $Q_a = 186.5$ kJ/mol and $n = 3.3$) and the gas constant $R (= 8.314$ J/kg/mol). The seed has constant viscosity of 10²¹ Pa · s. For models with a Cosserat continuum, the bending shear modulus and shear viscosity were kept equal to that used for the evaluation of deviatoric stress ($G_c = G$ and $\eta_c = \eta$). The remaining material parameters are listed in Table 1.

### 5. A Single Shear Band Development Model in 2D

#### 5.1. Verification of Cosserat Model Implementation Using FEM Simulations

The implementation of a Cosserat continuum in an existing code implies some fundamental changes (see Appendices A and B). We here verify our FDM implementation by benchmarking against FEM-based simulations.
<table><caption>**Table 1**

List of Material Parameters Employed in This Study</caption><thead><tr><th></th><th>G (Pa)</th><th>K (Pa)</th><th>C (Pa)</th><th>h (Pa)</th><th>φ (°)</th><th>ψ (°)</th><th>ε<sup>BG</sup> (s<sup>-1</sup>)</th><th>Δt (s)</th><th>R (m)</th></tr></thead><tbody><tr><td>Single shear band</td><td>1 × 10<sup>10</sup></td><td>2 × 10<sup>10</sup></td><td>1.75 × 10<sup>7</sup></td><td>-2 × 10<sup>8</sup></td><td>30</td><td>0-3-6</td><td>2 × 10<sup>-13</sup></td><td>8 × 10<sup>8</sup></td><td>500</td></tr><tr><td>Random shear bands</td><td>1 × 10<sup>10</sup></td><td>2 × 10<sup>10</sup></td><td>1.00 × 10<sup>8</sup></td><td>-2 × 10<sup>8</sup></td><td>30</td><td>0</td><td>2 × 10<sup>-13</sup></td><td>8 × 10<sup>8</sup></td><td>-</td></tr><tr><td>Crustal scale</td><td>1 × 10<sup>10</sup></td><td>2 × 10<sup>10</sup></td><td>1.75 × 10<sup>7</sup></td><td>-2 × 10<sup>8</sup></td><td>30</td><td>0</td><td>10<sup>-15</sup></td><td>2.5 × 10<sup>10</sup></td><td>2 × 10<sup>3</sup></td></tr></tbody></table>

The reference FEM solution has been generated using a standard displacement-rotation framework (Hageman et al., 2021). The spatial discretization was performed with cubic B-splines (Hughes et al., 2005), using 94 × 62 elements (element size 150 × 168 m). This resolution corresponds to that used for the FDM calculations in this section. We used a backward Euler scheme for the temporal discretization. One difference of note between the FDM implementation presented in this paper, and the reference FEM implementation is that the finite element solution required significantly smaller time increments to localize accurately: We have used 40 versus 400 time steps for the FDM and the reference FEM scheme, respectively.

The results depicted in Figure 1 were computed using the single shear band model configuration. We report the temporal evolution of the average of the second deviatoric stress invariant. The results were obtained for three different length scales ranging from 100 to 300 m, with and without cohesion softening. In all cases, we obtained good agreement between FDM and FEM results. The only noticeable mismatch was obtained for a length scale of 100 m including softening. Under these conditions, the shear band width reached a value which was too close to the spatial resolution of either (FEM or FDM) model.

## 5.2. Shear Banding With Viscoplasticity, Gradient Plasticity and in Cosserat Continuum

We investigated the single shear band model configuration using $\eta^{\text{vp}} = 10^{18}$ Pa·s for viscoplastic regularization, $K_g = 5.10^{12}$ Pa·m$^2$ for gradient plasticity and $l_c = 80$ m for the Cosserat continuum. The simulations have been run for 40 time steps up to a final time of ≈1.1 Kyr. The plastic strain maps exhibit a finite width band of increased strain localization for each regularization approach (Figures 2a–2c). For the parameters considered, the widths of the bands are comparable. In comparison with gradient and Cosserat implementations, the strain gradient appears sharper in the case of viscoplasticity. By analyzing the spatial distribution of $\tau_{II}$ (Figures 2d–2f), we observe that viscoplastic regularization exhibits small-scale shear band splitting. It should be noted the gradient implementation also exhibits a shear band splitting (Figure 2e), which is not observed in the case of a Cosserat regularization (Figure 2f). The Cosserat model clearly exhibits high stress regions on each side of the shear band (Figure 2f).

We have monitored the time evolution of the spatially averaged second invariant of deviatoric stress ($\bar{\tau}_{II} = 1/S \int \tau_{II} dS$ where S is the model's surface). Each regularization approach exhibits a different transient behavior (Figure 3). However, all models reach the same peak stress level (≈140 MPa) and a similar residual stress level (≈125 MPa). The viscoplastic model exhibits a delayed softening which reflects the viscous relaxation

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Verification of Cosserat model implementation using 2D shear banding simulations. The results of the considered finite difference method approach (crosses) were compared against finite element method results (solid lines). Various Cosserat length scales were considered (a) 100 m, (b) 200 m, and (c) 300 m. Models were computed with (red) and without cohesion softening (blue).</figcaption></figure>
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Shear banding pattern obtained with the various regularizations after 40 time steps ($t = 1.1$ Kyr). The upper row correspond to accumulated plastic strain and the lower row depicts the second invariant of the deviatoric stress tensor. Viscoplasticity (A, D, $\eta^{\text{VP}} = 10^{18}$ Pa·s), gradient plasticity (B, E, $K_g = 5 \times 10^{12}$ Pa·m$^2$) and the use of a Cosserat continuum (C, F, $l_c = 80$ m).</figcaption></figure>

of the overstress. The gradient plasticity exhibits the steepest softening. Each implementation shows convergence of the stress-time curve upon mesh refinement. The Cosserat implementation clearly outperforms the other approaches since almost no difference is noticeable between the results of the lowest and the highest tested resolution.

Profiles of the plastic strain ($\epsilon^p$, Figure 4a), the stress invariant ($\tau_{II}$, Figure 4b) and the pressure ($P$, Figure 4c) were probed across the shear bands. As expected from the 2D maps, viscoplasticity provides a plug flow profile for the plastic strain while the profiles resulting from the gradient and Cosserat models are more gradual (Figure 4a). All regularization approaches show a decreased value of the stress and the pressure inside the shear bands (Figures 4b

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Stress-time curves obtained using the three different regularization approaches. The average of the second invariant of the deviatoric stress tensor is reported ($\bar{\tau}_{II}$). For each regularization approach, three grid resolutions were employed (46 × 30, 94 × 62, 190 × 126 cells).</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Profiles orthogonal to the shear band (coordinate *x*') for (a) the plastic strain, (b) the second invariant of the deviatoric stress, and (c) the pressure.</figcaption></figure>

and 4c). For viscoplasticity and gradient plasticity, the development of multiples are noticeable on both stress and pressure profiles (Figures 4b and 4c). For Cosserat, the high stress rims observed in the 2D stress maps (Figure 2f) are now clearly observable along the stress profile (Figure 4b).

In order to verify that the latter effect is not an artifact of mesh resolution, we have run a grid convergence test. The resolution was increased from 46 × 30 up to 382 × 254 cells. The results clearly show convergence upon mesh refinement (Figure 5). The stress ($\tau_{II}$) in the rims saturates at $\approx$140 MPa with increasing resolution while the center of the shear band and the far-field reach values of 117 and 125 MPa, respectively. Moreover, independently computed FEM solutions (Figure 1) exhibited the same behavior.

Further differences between regularization schemes arise when considering the effect of dilatancy. Using the same model configuration, we have studied the effect of increasing the dilation angle (Figure 6). For all regularization,

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Convergence test for the Cosserat model. The profile of deviatoric stress orthogonal to the shear band (coordinate *x*') was measured for various mesh resolutions.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Shear banding patterns obtained using the different regularization approaches (VP: viscoplasticity with $\eta^{\text{vp}} = 10^{18}$ Pa·s, Grad.: gradient with $K_g$, Coss.: Cosserat with $l_c = 80$ m) and for different values of their characteristic parameter. A resolution of 190 × 126 cells was employed (d–f). The color map corresponds the value of accumulated strain.</figcaption></figure>

non-zero dilatancy tends to smear out strain localization. The intensity of plastic strain thus decreases within shear bands. Nevertheless, gradient and Cosserat models manage to maintained strain localization within shear bands (Figures 6f and 6i). In the case of viscoplastic regularization, plastic strain also occurs outside the shear bands. With increasing dilation angle, the plastic strain gradient reduces and localization an thus less pronounced (Figure 6c).

We studied the width of shear bands for each regularization approach as a function of the main regularization parame-
ter. Widths (W) were extracted from the plastic strain profiles and correspond to the bandwidth of the best Gaussian fit.
For all approaches, shear band widths of several 100 of meters were obtained; they are thus numerically resolved. For
viscoplasticity, the characteristic width is proportional to the inverse of the Kelvin element viscosity ($\eta^{\text{vp}}$, Figure 7a).
For gradient plasticity, shear band width scales with the square root of the gradient parameter ($K_g$, Figure 7b). For the
Cosserat model, the width depends linearly on the length scale $l_c$ (Figure 7c). This is another benefit of the Cosserat
model in comparison to the other approaches as it allows for a straightforward control of shear band widths.

<figure><img src="image_9.png" /><figcaption>**Figure 7.** Shear band width measured for the three regularization approaches. (a) Viscoplasticity, (b) gradient plasticity, and (c) Cosserat continuum.</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**Figure 8.** (a) Stress-strain curves obtained using the three different regularization approaches. Shear banding pattern obtained with viscoplasticity (b), gradient regularization (c) and Cosserat (d). The color map corresponds to the magnitude of accumulated strain.</figcaption></figure>

## 6. Shear Band Networks Arising From Initial Random Cohesion

We carried out simulations with a random initial cohesion for a viscoplasticity, gradient plasticity and Cosserat continuum.

We report the evolution of averaged stress in Figure 8a. Similar to the single shear band case, gradient plasticity provides the most abrupt stress drop while viscoplasticity provides the most gradual stress drop. For the considered set of parameters, the gradient and Cosserat models result in approximately the same residual load, ≈140 MPa, while the viscoplastic model reaches ≈155 MPa, which reflects the overstress. For each regularization approach, the shear banding pattern consists of one main thick shear band that reflects on the four model boundaries (Figure 8b). The location of reflection points differs in all cases. Gradient plasticity yields the narrowest band, while the viscoplasticity and Cosserat models are characterized by a similar width. In the case of viscoplasticity, the generation of multiple bands can clearly be observed, while gradient and Cosserat models remain localized in a single band.

The shear banding patterns obtained with the three approaches are remarkably stable upon mesh refinement (Figure 9). The viscoplastic model is the most affected by a variation of the model's resolution, however, this only impacts the details of the shear band multiples.

We additionally explored the role of the governing parameters for each regularization approach (Figure 10). As in the single shear band case (e.g., Figure 7), the shear band broadens with an increase of either $\eta^{\text{vp}}$, $K_g$, or $l_c$, while the intensity of the plastic strain decreases with increasing shear band width. For viscoplasticity we observe more prominent multiples for an increasing Kelvin-element viscosity. Minor multiples were also obtained with of gradient plasticity, but only for the lowest value of $K_g$. The Cosserat models are characterized by robust shear banding pattern, which is merely affected by a change of resolution or governing parameter and which is devoid of any multiple shear banding.

## 7. Crustal Scale

We have carried out simulations of shear banding at the crustal scale using the three regularization techniques. The results are depicted in Figure 11, which depicts maps of second strain rate invariant after 1.07 Myr. The models undergo horizontal compression and include a brittle-ductile transition as well as a free surface. We focused on shear banding that occurs during visco-elasto-plastic stress build up. Therefore, the effect of finite strain was not included and the model geometries did not evolve with time.

As for the previous test (Figure 10), we can clearly observe the effect of the regularization. For each approach, strain localization loses intensity with increasing the regularization parameter. In particular, we observe that main shear bands broaden and that surrounding shear band multiples smear out. Shear banding
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Comparison of shear banding patterns computed with the different regularization approaches and for different numerical resolutions. Low resolution (LR, 94 × 62 cells (a–c)). Medium resolution (MR, 190 × 126 cells (d–f)). High resolution (HR, 382 × 254 cells (g–i)). The color map corresponds to the value of accumulated strain.</figcaption></figure>

patterns can even exhibit striking similarities (e.g., Figures 10c and 10i). Again, multiples are much more present using viscoplastic regularization which exhibits high frequency shear bands at low values of η^VP (Figure 10a) and which resembles patterns obtained using non regularized frictional plasticity (i.e., local and rate-independent).

<figure><img src="image_12.png" /><figcaption>**Figure 10.** Shear banding patterns obtained using the different regularization approaches (VP: viscoplasticity, Grad.: gradient, Coss.: Cosserat) and for different values of their characteristic parameter. A resolution of 190 × 126 cells was employed (d–f). The color map corresponds to the value of accumulated strain.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Crustal shear banding modeled using the different regularization approaches (VP: viscoplasticity, Grad.: gradient, Coss.: Cosserat) and for different values of their characteristic parameter. A resolution of 510 × 382 cells was employed (d–f). The color map corresponds the value of the second strain rate invariant. Results are depicted after 1,350 time steps (≈1.07 Myr).</figcaption></figure>

## 8. Combining Benefits of Temporal and Spatial Regularization

In some cases, it may be beneficial to combine the effects of several regularization schemes. In particular when both effects of temporal and spatial regularization are desired. Such combination is easily achievable in the context of the employed PT solving strategy. We have investigated the combination of gradient plasticity and viscoplasticity as well as Cosserat continuum with viscoplasticity, using the models based on the single shear band configuration. The stress-time evolution curves indicated the effect of viscoplastic overstress on the gradient and Cosserat model results (Figure 12a). The latter is minimized by setting low values of regularization viscosity ($\eta^{\text{VP}} < 5.10^{17}$ Pa·s). As observed for direct iterative schemes (e.g., Duretz et al., 2019), the inclusion of rate dependence also facilitates convergence to equilibrium. This is particularly remarkable for gradient plasticity (Figure 12b) where even small value of $\eta^{\text{VP}}$ can lead to 50% reduction in the number of iterations. The effect of viscoplasticity is also noticeable for the Cosserat model but the improvement is less than 25% (Figure 12c).

## 9. Discussion: The Benefits and Limitations of Regularization Approaches

For the configurations considered, each regularization approach delivers satisfactory results. They allow for convergence regarding mesh refinement, satisfaction of equilibrium, and result in comparable shear band widths. Nevertheless, each approach influences the solution in its own way.

The inclusion of viscoplasticity induces overstress and large values of the viscosity parameter, $\eta^{\text{VP}}$, can lead to an overshoot of the material strength. Nevertheless, for moderate overstress, fairly accurate strength estimates can be achieved while keeping the regularization benefits of viscoplasticity (Duretz et al., 2020). A power-law model can also be used to limit the effects of growing overstress with increasing strain rate (Duretz et al., 2021). Another drawback of viscoplasticity is the occurrence of multiple shear bands and a possible delocalization of the strain. The latter is particularly well expressed when non-zero dilatancy is taken into account. This effect was vastly reduced with gradient plasticity and was not observed with models based on a Cosserat continuum. In contrast with viscoplasticity and gradient plasticity, in the Cosserat continuum, shear bands exhibit rims which of high differential stress, which likely contributes to the focusing of the strain inside the band. A comparison with stress distributions around natural shear/fault zones could help to validate this. In general the Cosserat model provides the most robust shear banding patterns and is the simplest way to control shear band width. At this point, we only
<figure><img src="image_14.png" /><figcaption>**Figure 12.** Combination of spatial and temporal regularizations. Panel a depicts stress-time evolution curves. Panels b and c show the pseudo-transient iteration count for each time step for gradient/viscoplasticity and Cosserat/viscoplasticity, respectively. On each panel the blue line corresponds to the results obtained solely with viscoplasticity and using the reference value $\eta^{\text{vp}} = 10^{18} \text{ Pa} \cdot \text{s}$. Similarly, gradient and Cosserat parameters were set to $K_g = 5.10^{12} \text{ Pa} \cdot \text{m}^2$ and $l_c = 80 \text{ m}$.</figcaption></figure>

intended to use Cosserat as a numerical regularization technique. We have made simplifying assumptions in the choice of the extra parameters involved in the Cosserat model. Notably, the definition of invariants involves extra material parameters. By assuming $h_1 = h_2 = h_3$, we have restricted our models to isotropic viscous and plastic flow. Accurate constitutive modeling using Cosserat medium can also be achieved (e.g., Stefanou et al., 2017). Such models often imply mechanical anisotropy as reflected by estimated values $h_1, h_2$, and $h_3$ for different materials.

# 10. Implications and Perspectives for Geodynamic Modeling

The considered regularization approaches all show attractive properties for geodynamic models which typically involve the development of fault zones, shear zones and plate boundaries. It should be noted that the considered regularizations do not alter the angle of shear bands arising from shear instabilities using Drucker-Prager rheology. Spatial regularizations (gradient and Cosserat continuum) provide the most robust shear banding patterns. They manage to maintain plastic deformation within the shear bands, exhibit excellent mesh convergence and simple control over shear band width.

Although we have used a Cosserat continuum as a regularization, with an appropriate description of the material structure, a Cosserat model can adequately reflect the mechanical behavior of a material (e.g., Stefanou et al., 2017). Moreover, Cosserat continua are general because, besides plasticity, they also introduce a length scale in elasticity (e.g., Lakes, 1995) or viscosity (Riihimäki, 1978). The latter could have interesting implications in geophysics such as in seismology (Abreu et al., 2017). Upon adequate upscaling of rocks microstructure, the use of Cosserat media could provide further insights in the patterns of interseismic deformation, which are typically studied using
classical elasticity and viscosity (Savage & Burford, 1973; Traoré et al., 2014). The role of couple-stresses was shown to be beneficial for the simulation of anisotropic multilayer folding (Muhlhaus et al., 2002). Additionally, Cosserat media may constitute an appealing alternative for modeling ductile strain localization using a constitutive model that contains no inherent length scale or length scales that are typically non-resolvable in a typical geodynamic simulation (e.g., grain size evolution, Rozel et al., (2011)).

## 11. Implications for Geodynamic Simulation Tools

Viscoplasticity is straightforward to implement and is already used in codes such as LaMEM (Moulas et al., 2022), M2Di (Duretz et al., 2020) and MDoodz (Duretz et al., 2021; Yamato et al., 2022). Nevertheless, a finer control on shear banding patterns is achievable using spatial regularizations such as gradient plasticity or Cosserat continuum. In particular, the use of a Cosserat continuum allows for a direct control of shear band widths (e.g., (Mühlhaus & Vardoulakis, 1987), see also Figure 7c). However the implementation of spatial regularization requires major changes in the structure and solving strategies of existing codes. With gradient regularization, the return mapping procedure involves the Laplacian of the plastic multiplier rate. This involves the solution of an additional coupled partial differential equation with boundary conditions defined at the boundary of the plastic domain. The plastic multiplier rate thus becomes an additional global DOF (in 2D or 3D). In the Cosserat continuum, the components of the micro-rotation rate vector are additional DOFs. This corresponds to one and three additional global DOF, and associated partial differential equation, in 2D and 3D, respectively.

The PT solving strategy enabled us to analyze the different regularization in a single framework. Nevertheless, the inclusion of spatial regularization techniques may also be easily achieved with other solving strategies that are designed to automatically discretize and solve coupled non-linear systems of equations (Davies et al., 2022; Wilson et al., 2017). Alternatively, it may also be possible to include effects of spatial regularization by decoupling the equations (operator split) and using staggered schemes such as usually applied in damage modeling (Miehe et al., 2010).

## 12. Conclusions

Frictional strain localization is an important feature of geodynamic simulations. However, due to the lack of inherent length and time scales, the numerical treatment of frictional plasticity in geodynamic model is problematic. For this reason, we have investigated the effects of three temporal and spatial regularization techniques: viscoplasticity, gradient plasticity and the use of a Cosserat continuum. The three techniques were implemented in a single code based on the accelerated PT method combined to a FDM discretization. To conform to geodynamic modeling standards, we expressed the Cosserat model for velocity-pressure formulations and benchmarked it using state-of-the-art FEM models. The three regularization strategies provide attractive properties for modeling strain localization in frictional material equipped (or not) with a material strain softening parametrization. All regularization strategies allow, to some extent, for controlling shear band widths, they deliver convergence upon mesh refinement and satisfaction of the force balance. All regularization techniques were successful at modeling crustal scale shear banding during visco-elasto-plastic build-up. In practice, viscoplasticity is straightforward to include in existing codes. Gradient and Cosserat implementations require extra degrees of freedom but both provide more robust shear banding patterns and more precise strength estimates. Besides regularization, the use of Cosserat-type media to upscale rock fabrics brings exciting perspectives for future geodynamic modeling.

## Appendix A: Splitting of the Shear Stress Components

In the Cosserat model, the total strain rate tensor is not symmetric as the shear components include effect of micro-rotation rates:

$$ 
\begin{aligned}
\dot{\epsilon}_{xy} &= \frac{\partial v_y}{\partial x} - \dot{\omega}_z \\
\dot{\epsilon}_{yx} &= \frac{\partial v_x}{\partial y} + \dot{\omega}_z
\end{aligned}
\qquad (A1) $$

Assuming a viscous rheology, the total shear stress component $\sigma_{xy} \neq \sigma_{yx}$ and is expressed as:

$$ 
\begin{aligned}
\sigma_{xy} &= \eta \dot{\epsilon}_{xy} + \eta_c \dot{\epsilon}_{xy} + \eta \dot{\epsilon}_{yx} - \eta_c \dot{\epsilon}_{yx} \\
\sigma_{yx} &= \eta \dot{\epsilon}_{yx} + \eta_c \dot{\epsilon}_{yx} + \eta \dot{\epsilon}_{xy} - \eta_c \dot{\epsilon}_{xy}
\end{aligned}
\qquad (A2) $$
Upon expansion, $\sigma_{xy}$ may be reformulated as:

$$ 
\begin{aligned}
\sigma_{xy} &= \eta \frac{\partial v_y}{\partial x} - \eta \dot{\omega}_z + \eta_c \frac{\partial v_y}{\partial x} - \eta_c \dot{\omega}_z + \eta \frac{\partial v_x}{\partial y} + \eta \dot{\omega}_z + \eta_c \frac{\partial v_x}{\partial y} - \eta_c \dot{\omega}_z \\
\sigma_{xy} &= \eta \left( \frac{\partial v_y}{\partial x} + \frac{\partial v_x}{\partial y} \right) + \eta_c \left( \frac{\partial v_y}{\partial x} - \frac{\partial v_x}{\partial y} \right) - 2\eta_c \dot{\omega}_z
\end{aligned}
\quad (A3) $$

or after using the symmetrized shear strain rate $\dot{\epsilon}_{xy} = \frac{1}{2} \left( \frac{\partial v_y}{\partial x} + \frac{\partial v_x}{\partial y} \right)$ and the spin $\dot{\omega}_z = \frac{1}{2} \left( \frac{\partial v_x}{\partial y} - \frac{\partial v_y}{\partial x} \right)$

$$ \sigma_{xy} = 2\eta\dot{\epsilon}_{xy} - 2\eta_c\dot{\omega}_z - 2\eta_c\dot{\omega}_z \quad (A4) $$

Using the identities $\tau_{xy} = 2\eta\dot{\epsilon}_{xy}$ and that $R_z = -2\eta_c(\dot{\omega}_z + \dot{\omega}_x)$, and carrying out a similar derivation for $\sigma_{yx}$, one obtains:

$$ \sigma_{xy} = \tau_{xy} + R_z \quad (A5) $$

$$ \sigma_{yx} = \tau_{xy} - R_z $$

These expressions are convenient to work with in geodynamic code since the deviatoric shear stress remains symmetric. The only addition is that one needs to keep track of the additional (antisymmetric) components of shear stress caused by the rotation ($R_z$) in order to calculate the total shear ($\sigma_{yx}$) and subsequently evaluate the linear and angular momentum balances.

## Appendix B: Cosserat Model for Velocity-Pressure Formulation

Traditionally implementations of Cosserat are expressed in a total stress-displacement formulation and based on elasto-plastic rheological model (e.g., de Borst, 1991). Here we express the Cosserat continuum model in terms of velocity-pressure formulation and based on a visco-elasto-plastic rheological model. We seek a solution the following system of coupled equations:

$$ 
\begin{aligned}
-\frac{\partial p}{\partial x} + \frac{\partial \tau_{xx}}{\partial x} + \frac{\partial \tau_{xy}}{\partial y} + \frac{\partial R_z}{\partial y} &= 0 \\
-\frac{\partial p}{\partial y} + \frac{\partial \tau_{yy}}{\partial y} + \frac{\partial \tau_{xy}}{\partial x} - \frac{\partial R_z}{\partial x} &= 0 \\
\frac{\partial m_{xz}}{\partial x} + \frac{\partial m_{yz}}{\partial y} + 2R_z &= 0 \\
\frac{\dot{p}}{K} + \frac{\partial v_x}{\partial x} + \frac{\partial v_y}{\partial y} &= 0
\end{aligned}
\quad (B1) $$

The solutions, or primitive variable, are: $v_x, v_y, \dot{w}_z$, and $p$. The balance angular momentum shows that the total stress tensor is non-symmetric and the shear stress components may be expressed as: $\sigma_{xy} = \tau_{xy} + R_z$ and $\sigma_{yx} = \tau_{xy} - R_z$ (see above).

The following kinematic relationships relate volumetric, deviatoric, and rotational strain rate to the primitive variables:

$$ 
\begin{aligned}
\nabla \cdot v &= \frac{\partial v_x}{\partial x} + \frac{\partial v_y}{\partial y} \\
\dot{\epsilon}'_{xx} &= -\frac{1}{3} \nabla \cdot v + \frac{\partial v_x}{\partial x} \\
\dot{\epsilon}'_{yy} &= -\frac{1}{3} \nabla \cdot v + \frac{\partial v_y}{\partial y} \\
\dot{\epsilon}'_{zz} &= -\frac{1}{3} \nabla \cdot v \\
\dot{\epsilon}'_{xy} &= \frac{1}{2} \left( \frac{\partial v_x}{\partial y} + \frac{\partial v_y}{\partial x} \right) \\
\dot{w}_z &= \frac{1}{2} \left( \frac{\partial v_x}{\partial y} - \frac{\partial v_y}{\partial x} \right) \\
\dot{k}_{xz} &= \frac{\partial \dot{\omega}_z}{\partial x} \\
\dot{k}_{yz} &= \frac{\partial \dot{\omega}_z}{\partial y}
\end{aligned}
\quad (B2) $$
where the notation $\dot{\epsilon}'_{ij}$ indicate the symmetrized deviatoric strain rate tensor components. Plane-strain conditions are assumed. The constitutive relationships are expressed in the continuous form in Section 2.2. In the following, we employ a semi-discrete form for practical purposes. The rates are discretized using backward Euler time integration. We assume small elastic strains, such that terms related to stress tensor advection, rotation and stretch are negligible. The updates of deviatoric stress, force-conjugate of micro-rotation rate and couple-stress may be expressed as:

$$ 
\begin{aligned}
\tau_{xx} &= 2\eta^{\text{ve}} \left( \dot{\epsilon}'_{xx} + \frac{\tau_{xx}^0}{2G\Delta t} \right) \\
\tau_{yy} &= 2\eta^{\text{ve}} \left( \dot{\epsilon}'_{yy} + \frac{\tau_{yy}^0}{2G\Delta t} \right) \\
\tau_{zz} &= 2\eta^{\text{ve}} \left( \dot{\epsilon}'_{zz} + \frac{\tau_{zz}^0}{2G\Delta t} \right) \\
\tau_{xy} &= 2\eta^{\text{ve}} \left( \dot{\epsilon}'_{xy} + \frac{\tau_{xy}^0}{2G\Delta t} \right) \\
R_z &= -2\eta_c^{\text{ve}} \left[ (\dot{w}_z + \dot{\omega}_z) - \frac{R_z^0}{2G_c\Delta t} \right] \\
m_{xz} &= l_c^2 2\eta^{\text{ve}} \left( \dot{k}_{xz} + \frac{m_{xz}^0}{2l_c^2 G\Delta t} \right) \\
m_{yz} &= l_c^2 2\eta^{\text{ve}} \left( \dot{k}_{yz} + \frac{m_{yz}^0}{2l_c^2 G\Delta t} \right)
\end{aligned} 
\quad (B3) $$

where $\eta^{\text{ve}} = \left(\frac{1}{\eta} + \frac{1}{G\Delta t}\right)^{-1}$ is the effective modulus for deviatoric visco-elasticity and $\eta_c^{\text{ve}} = \left(\frac{1}{\eta_c} + \frac{1}{G_c\Delta t}\right)^{-1}$ is the effective modulus for Cosserat rotation visco-elasticity. The rotation component $R_z$ is anti-symmetric and the deviatoric shear stress is symmetric $\tau_{xy} = \tau_{yx}$. The superscript 0 indicates quantities from the previous time step and result from the backward-Euler integration. The second deviatoric stress invariant takes the form of:

$$ \tau_{II} = \sqrt{h_1 \tau_{xx}^2 + h_1 \tau_{yy}^2 + h_1 \tau_{zz}^2 + 2h_2 \tau_{xy}^2 + 2h_2 R_z^2 + h_3 \left(\frac{m_{xz}}{l_c}\right)^2 + h_3 \left(\frac{m_{yz}}{l_c}\right)^2} \quad (B4) $$

For the Drucker-Prager model, we define the plastic flow potential, $Q = \tau_{II} - p \sin \psi$. At trial state, we define $p^{\text{trial}} = p$. The partial derivatives of the plastic flow potential are:

$$ 
\begin{aligned}
\frac{\partial Q}{\partial \tau_{xx}} &= h_1 \frac{\tau_{xx}}{\tau_{II}} \\
\frac{\partial Q}{\partial \tau_{yy}} &= h_1 \frac{\tau_{yy}}{\tau_{II}} \\
\frac{\partial Q}{\partial \tau_{zz}} &= h_1 \frac{\tau_{zz}}{\tau_{II}} \\
\frac{\partial Q}{\partial \tau_{xy}} &= 2h_2 \frac{\tau_{xy}}{\tau_{II}} \\
\frac{\partial Q}{\partial R_z} &= 2h_2 \frac{R_z}{\tau_{II}} \\
\frac{\partial Q}{\partial m_{xz}} &= h_3 \frac{m_{xz}}{\tau_{II} l_c^2} \\
\frac{\partial Q}{\partial m_{yz}} &= h_3 \frac{m_{yz}}{\tau_{II} l_c^2} \\
\frac{\partial Q}{\partial p} &= -\sin \psi
\end{aligned} 
\quad (B5) $$

Finally, the corrected deviatoric stresses read:
$$
\begin{align*}
\tau_{xx} &= 2\eta^{\text{ve}} \left( \dot{\epsilon}_{xx} + \frac{\tau_{xx}^0}{2G\Delta t} - \dot{\lambda} \frac{\partial Q}{\partial \tau_{xx}} \right) \\
\tau_{yy} &= 2\eta^{\text{ve}} \left( \dot{\epsilon}_{yy} + \frac{\tau_{yy}^0}{2G\Delta t} - \dot{\lambda} \frac{\partial Q}{\partial \tau_{yy}} \right) \\
\tau_{zz} &= 2\eta^{\text{ve}} \left( \dot{\epsilon}_{zz} + \frac{\tau_{zz}^0}{2G\Delta t} - \dot{\lambda} \frac{\partial Q}{\partial \tau_{zz}} \right) \\
\tau_{xy} &= 2\eta^{\text{ve}} \left( \dot{\epsilon}_{xy} + \frac{\tau_{xy}^0}{2G\Delta t} - \dot{\lambda} \frac{\partial Q}{\partial \tau_{xy}} \right) \\
R_z &= -2\eta_c^{\text{ve}} \left( (\dot{w}_{xy} + \dot{\omega}_z) - \frac{R_z^0}{2G_c\Delta t} - \dot{\lambda} \frac{\partial Q}{\partial R_z} \right) \\
m_{xz} &= 2\eta^{\text{ve}} l_c^2 \left( \dot{k}_{xz} + \frac{m_{xz}^0}{2G\Delta t l_c^2} - \frac{\dot{\lambda}}{l_c} \frac{\partial Q}{\partial m_{xz}} \right) \\
m_{yz} &= 2\eta^{\text{ve}} l_c^2 \left( \dot{k}_{yz} + \frac{m_{yz}^0}{2G\Delta t l_c^2} - \frac{\dot{\lambda}}{l_c} \frac{\partial Q}{\partial m_{yz}} \right) \\
p &= p^{\text{trial}} - K\Delta t \dot{\lambda} \frac{\partial Q}{\partial p} \tag{B6}
\end{align*}
$$

The plastic multiplier rate, $\dot{\lambda}$, is expressed as

$$ \dot{\lambda} = \frac{F^{\text{trial}}}{\eta^{\text{ve}} + \eta^{\text{vp}} + K \Delta t \sin \phi \sin \psi + \Delta t \sqrt{\frac{z}{3}} h \cos \phi} \quad (B7) $$

where $F^{\text{trial}}$ is the yield function evaluated for trial deviatoric stress invariant and trial pressure. Furthermore, if one assumes that $\eta_c^{\text{ve}} = \eta^{\text{ve}}$, the deviatoric stress invariant becomes proportionally to the effective strain rate invariant (see e.g., Kaus et al., 2016; Moresi et al., 2003) and an effective viscosity can be expressed (e.g., Duretz et al., 2021).

## Data Availability Statement

All the scripts are available in the repository (Räss & Duretz, 2023), the initial cohesion field is available in this repository (Duretz, 2023). The scripts are written in the open source Julia language and the scripts' latest version is accessible at https://github.com/PTsolvers/PlasticityRegularisations_G3.

## References

Abreu, R., Kamm, J., & Reiß, A.-S. (2017). Micropolar modelling of rotational waves in seismology. *Geophysical Journal International*, 210(2), 1021–1046. https://doi.org/10.1093/gji/ggx211

Bažant, Z. P., & Lin, F.-B. (1988). Non-local yield limit degradation. *International Journal for Numerical Methods in Engineering*, 26(8), 1805–1823. https://doi.org/10.1002/nme.1620260809

Bezanson, J., Edelman, A., Karpinski, S., & Shah, V. B. (2017). Julia: A fresh approach to numerical computing. *SIAM Review*, 59(1), 65–98. https://doi.org/10.1137/141000671

Brantut, N., Stefanou, I., & Sulem, J. (2017). Dehydration-induced instabilities at intermediate depths in subduction zones. *Journal of Geophysical Research: Solid Earth*, 122(8), 6087–6107. https://doi.org/10.1002/2017JB014357

Buiters, S. J., Schreurs, G., Alberts, M., Gerya, T. V., Kaus, B., Landry, W., et al. (2016). Benchmarking numerical models of brittle thrust wedges. *Journal of Structural Geology*, 92, 140–177. https://doi.org/10.1016/j.jsg.2016.03.003

Byerlee, J. (1978). Friction of rocks. *Pure and Applied Geophysics*, 116(4–5), 615–626. https://doi.org/10.1007/bf00876528

Davies, D. R., Kramer, S. C., Ghelichkhan, S., & Gibson, A. (2022). Towards automatic finite-element methods for geodynamics via firedrake. *Geoscientific Model Development*, 15(13), 5127–5166. https://doi.org/10.5194/gmd-15-5127-2022

de Borst, R. (1991). Simulation of strain localization: A reappraisal of the Cosserat continuum. *Engineering Computations*, 8(4), 317–332. https://doi.org/10.1108/eb023842

de Borst, R., & Duretz, T. (2020). On viscoplastic regularisation in strain-softening rocks and soils. *International Journal for Numerical and Analytical Methods in Geomechanics*, 44(6), 890–903. https://doi.org/10.1002/nag.3046

de Borst, R., & Mühlhaus, H. B. (1992). Gradient-dependent plasticity: Formulation and algorithmic aspects. *International Journal for Numerical Methods in Engineering*, 35(3), 521–539. https://doi.org/10.1002/nme.1620350307

de Borst, R., Sluys, L. J., Mühlhaus, H.-B., & Pamin, J. (1993). Fundamental issues in finite element analysis of localization of deformation. *Engineering Computations*, 10(2), 99–122. https://doi.org/10.1108/eb023897

Duretz, T. (2023). Random cohesion field [Dataset]. Zenodo. https://doi.org/10.5281/zenodo.7692569

Duretz, T., de Borst, R., & Le Pourhiet, L. (2019). On finite thickness of shear bands in frictional viscoplasticity, and implications for lithosphere dynamics. *Geochemistry, Geophysics, Geosystems*, 20(11), 5598–5616. https://doi.org/10.1029/2019gc008531

## Acknowledgments

The authors thank Daniel Kiss and an anonymous reviewer for their constructive suggestions. TD acknowledges Rennes Métropole and the French programme NEEDS for their financial supports to the project LR acknowledges the Swiss Geocomputing Centre, University of Lausanne, and the Laboratory of Hydraulics, Hydrology and Glaciology (VAW), ETH Zurich, for providing computing resources on the Octopus and Superzack GPU-powered supercomputer and server, respectively. Open Access funding enabled and organized by Projekt DEAL.
Duretz, T., de Borst, R., & Yamato, P. (2021). Modeling lithospheric deformation using a compressible visco-elasto-viscoplastic rheology and the effective viscosity approach. *Geochemistry, Geophysics, Geosystems*, 22(8), e2021GC009675. https://doi.org/10.1029/2021GC009675

Duretz, T., de Borst, R., Yamato, P., & Le Pourhiet, L. (2020). Toward robust and predictive geodynamic modeling: The way forward in frictional plasticity. *Geophysical Research Letters*, 47(5), e2019GL086027. https://doi.org/10.1029/2019GL086027

Duretz, T., Räss, L., Podladchikov, Y., & Schmalholz, S. (2018). Resolving thermomechanical coupling in two and three dimensions: Spontaneous strain localization owing to shear heating. *Geophysical Journal International*, 216(1), 365–379. https://doi.org/10.1093/gji/ggy434

Duretz, T., Souche, A., de Borst, R., & Le Pourhiet, L. (2018). The benefits of using a consistent tangent operator for viscoelastoplastic computations in geodynamics. *Geochemistry, Geophysics, Geosystems*, 19(12), 4904–4024. https://doi.org/10.1029/2018gc007877

Gerya, T. (2013). Initiation of transform faults at rifted continental margins: 3D petrological-thermomechanical modeling and comparison to the woodlark basin. *Petrology*, 21(6), 550–560. https://doi.org/10.1134/S0869591113060039

Hageman, T., Sabet, S. A., & de Borst, R. (2021). Convergence in non-associated plasticity and fracture propagation for standard, rate-dependent, and Cosserat continua. *International Journal for Numerical Methods in Engineering*, 122(3), 777–795. https://doi.org/10.1002/nme.6561

Hansen, F., & Carter, N. (1983). Semibrittle creep of dry and wet westerly granite at 1000 MPa.

Heeres, O. M., Suiker, A. S. J., & de Borst, R. (2002). A comparison between the Perzyna viscoplastic model and the consistency viscoplastic model. *European Journal of Mechanics—A: Solids*, 21, 1–12. https://doi.org/10.1016/s0997-7538(01)01188-3

Hughes, T. J. R., Cottrell, J. A., & Bazilevs, Y. (2005). Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement. *Computer Methods in Applied Mechanics and Engineering*, 194(39–41), 4135–4195. https://doi.org/10.1016/j.CMA.2004.10.008

Jacquey, A. B., & Cacace, M. (2020). Multiphysics modeling of a brittle-ductile lithosphere: 1. Explicit visco-elasto-plastic formulation and its numerical implementation. *Journal of Geophysical Research: Solid Earth*, 125(1), e2019JB018474. https://doi.org/10.1029/2019jb018474

Kaus, B. J. P., Popov, A. A., Baumann, T. S., Püsök, A. E., Bauville, A., Fernandez, N., & Collignon, M. (2016). Forward and inverse modelling of lithospheric deformation on geological timescales. In K. Binder, M. Müller, A. Kremer, & A. Schnurpfeil (Eds.), *NIC symposium 2016* (Vol. 48, pp. 299–307). Forschungszentrum Jülich.

Lakes, R. S. (1995). Experimental methods for study of Cosserat elastic solids and other generalized elastic continua.

Leuthold, J., Gerolymatou, E., & Triantafyllidis, T. (2021). Effect of compaction banding on the hydraulic properties of porous rock—Part II: Constitutive description and numerical simulations. *Rock Mechanics and Rock Engineering*, 54(6), 2685–2696. https://doi.org/10.1007/s00603-021-02430-1

Miehe, C., Hofacker, M., & Welschinger, F. (2010). A phase field model for rate-independent crack propagation: Robust algorithmic implementation based on operator splits. *Computer Methods in Applied Mechanics and Engineering*, 199(45), 2765–2778. https://doi.org/10.1016/j.cma.2010.04.011

Moresi, L., Dufour, F., & Muhlhaus, H. (2003). A Lagrangian integration point finite element method for large deformation modeling of viscoelastic geomaterials. *Journal of Computational Physics*, 184(2), 476–497. https://doi.org/10.1016/s0021-9991(02)00031-1

Moulas, E., Kaus, B., & Jamtveit, B. (2022). Dynamic pressure variations in the lower crust caused by localized fluid-induced weakening. *Communications Earth & Environment*, 3(1), 1234567890. https://doi.org/10.1038/s43247-022-00478-7

Muhlhaus, H.-B., Dufour, F., Moresi, L., & Hobbs, B. E. (2002). A director theory for visco-elastic folding instabilities in multilayered rock. *International Journal of Solids and Structures*, 39(13–14), 3675–3691. https://doi.org/10.1016/s0020-7683(02)00175-0

Mühlhaus, H. B., & Vardoulakis, I. (1987). The thickness of shear bands in granular materials. *Géotechnique*, 37(3), 271–283. https://doi.org/10.1680/geot.1987.37.3.271

Naliboff, J. B., Buiter, S. J. H., Péron-Pinvidic, G., Osmundsen, P. T., & Tetreault, J. (2017). Complex fault interaction controls continental rifting. *Nature Communications*, 8(1), 1179. https://doi.org/10.1038/s41467-017-00904-x

Omlin, S., & Räss, L. (2021a). ImplicitGlobalGrid.jl.

Omlin, S., & Räss, L. (2021b). ParallelStencil.jl.

Pranger, C., Sanan, P., May, D., Pourhiet, L., & Gabriel, A. (2022). Rate and state friction as a spatially regularized transient viscous flow law. *Journal of Geophysical Research: Solid Earth*, 127(6), e2021JB023511. https://doi.org/10.1029/2021JB023511

Räss, L., & Duretz, T. (2023). PTsolvers/PlasticityRegularisations_G3: Initial release (v1.0.0). Zenodo. https://doi.org/10.5281/zenodo.7971379

Räss, L., Duretz, T., & Podladchikov, Y. Y. (2019). Resolving hydromechanical coupling in two and three dimensions: Spontaneous channelling of porous fluids owing to decompaction weakening. *Geophysical Journal International*, 218(3), 1591–1616. https://doi.org/10.1093/gji/ggz239

Räss, L., Duretz, T., Podladchikov, Y. Y., & Schmalholz, S. M. (2017). M2Di: Concise and efficient MATLAB 2-D Stokes solvers using the finite difference method. *Geochemistry, Geophysics, Geosystems*, 18(2), 755–768. https://doi.org/10.1002/2016gc006727

Räss, L., Utkin, I., Duretz, T., Omlin, S., & Podladchikov, Y. Y. (2022). Assessing the robustness and scalability of the accelerated pseudo-transient method. *Geoscientific Model Development*, 15(14), 5757–5786. https://doi.org/10.5194/gmd-15-5757-2022

Riihimäki, H. (1978). The Cosserat continuum and a viscous fluid. *Physics Letters A*, 68(2), 189–190. https://doi.org/10.1016/0375-9601(78)90799-5

Rozel, A., Ricard, Y., & Bercovici, D. (2011). A thermodynamically self-consistent damage equation for grain size evolution during dynamic recrystallization. *Geophysical Journal International*, 184(2), 719–728. https://doi.org/10.1111/j.1365-246X.2010.04875.x

Rudnicki, J. W., & Rice, J. R. (1975). Conditions for the localization of deformation in pressure sensitive dilatant materials. *Journal of the Mechanics and Physics of Solids*, 23(6), 371–394. https://doi.org/10.1016/0022-5096(75)90001-0

Sabet, S. A., & de Borst, R. (2019). Structural softening, mesh dependence, and regularisation in non-associated plastic flow. *International Journal for Numerical and Analytical Methods in Geomechanics*, 43(13), 2170–2183. https://doi.org/10.1002/nag.2973

Savage, J. C., & Burford, R. O. (1973). Geodetic determination of relative plate motion in central California. *Journal of Geophysical Research*, 78(5), 832–845. https://doi.org/10.1029/JB078i005p00832

Schmalholz, S. M., Moulas, E., Plümper, O., Myasnikov, A. V., & Podladchikov, Y. Y. (2020). 2D hydro-mechanical-chemical modeling of (de) hydration reactions in deforming heterogeneous rock: The periclase-brucite model reaction. *Geochemistry*, 21(11), 2020GC009351. https://doi.org/10.1029/2020gc009351

Spiegelman, M., May, D. A., & Wilson, C. R. (2016). On the solvability of incompressible Stokes with viscoplastic rheologies in geodynamics. *Geochemistry, Geophysics, Geosystems*, 17(6), 2213–2238. https://doi.org/10.1002/2015gc006228

Stathas, A., & Stefanou, I. (2021). The role of viscous regularization in dynamical problems, strain localization and mesh dependency. *Computer Methods in Applied Mechanics and Engineering*, 388, 114185. https://doi.org/10.1016/j.cma.2021.114185

Stefanou, I., Sulem, J., & Rattez, H. (2017). Cosserat approach to localization in geomaterials (pp. 1–25). Springer International Publishing. https://doi.org/10.1007/978-3-319-22977-5_10-1

Traoré, N., Le Pourhiet, L., Frelat, J., Rolandone, F., & Meyer, B. (2014). Does interseismic strain localization near strike-slip faults result from boundary conditions or rheological structure? *Geophysical Journal International*, 197(1), 50–62. https://doi.org/10.1093/gji/ggu011

Wang, H. (2019). Viscous and second gradient regularization techniques for the description of the behavior of geomaterials Theses. École centrale de Nantes.
Wang, W. M., Sluys, L. J., & de Borst, R. (1997). Viscoplasticity for instabilities due to strain softening and strain-rate softening. *International Journal for Numerical Methods in Engineering, 40*(20), 3839–3864. https://doi.org/10.1002/(sici)1097-0207(19971030)40:20<3839::aid-nme245>3.0.co;2-6

Wilson, C. R., Spiegelman, M., & van Keken, P. E. (2017). TerraFERMA: The transparent finite element rapid model assembler for multiphysics problems in Earth sciences. *Geochemistry, Geophysics, Geosystems, 18*(2), 769–810. https://doi.org/10.1002/2016GC006702

Yamato, P., Duretz, T., Baïsset, M., & Luisier, C. (2022). Reaction-induced volume change triggers brittle failure at eclogite facies conditions. *Earth and Planetary Science Letters, 584*, 117520. https://doi.org/10.1016/j.epsl.2022.117520

Zhao, X., & Cai, M. (2010). Influence of plastic shear strain and confinement-dependent rock dilation on rock failure and displacement near an excavation boundary. *International Journal of Rock Mechanics and Mining Sciences, 47*(5), 723–738. https://doi.org/10.1016/j.ijrmms.2010.04.003