A dilatant visco-elasto-viscoplasticity model with globally continuous tensile cap: stable two-field mixed formulation

Rocks break if shear stresses exceed their strength. It is therefore important for typical geoscientific applications to take shear failure mechanism and the subsequent development of mode-II shear bands or faults into account. Many existing codes incorporate non-associated Drucker-Prager or Mohr-Coulomb plasticity models to simulate this behavior. Yet, when effective mean stress becomes extensional, for example when fluid pressure becomes large, the dominant failure mode changes to a mode-I (opening) mode, which initiates plastic volumetric deformation. It is rather difficult to represent both failure modes in numerical models in a self-consistent manner, while also accounting for the nonlinear visco-elastic host rock rheology, which varies from being nearly incompressible in the mantle to being compressible in surface-near regions. Here, we present a simple plasticity model that is designed to overcome these difficulties. We employ a combination of a linearized Drucker-Prager shear failure envelope with a circular tensile cap function in way that ensures continuity and smoothness of both yield surface and flow potential in the entire stress space. A Perzyna-type viscoplastic regularization ensures that the resulting localization zones are mesh-insensitive. To deal with the near incompressibility condition, a mixed two-field finite element formulation is employed. The local nonlinear iterations at the integration-point level are used to determine the stress increments. The global Newton-Raphson iterations are applied to solve the discretized momentum and continuity residual equations. The presented plasticity model is implemented in an open-source 2D unstructured finite element code GeoTech2D. The results of several typical test cases that range from crustal scale deformation to the propagation of fluid-induced tensile failure zones demonstrate rapid convergence. The robustness of the solution scheme is enhanced by the adaptive time stepping algorithm.

1 Introduction

Plastic deformation in brittle rocks manifests itself in two ways: as mode-II shear faults, or as mode-I failure zones, which results in an opening of the rock in a crack-like manner and therefore produces volumetric strains. Having the ability to simulate both faulting modes numerically in a self-consistent manner is important to simulate near-surface brittle deformation, or simulate cases where the fluid pressure is high, such as during the initiation of hydraulic fractures or magma-filled dykes.

Quasi-static long-term flow of viscous materials, such as hot rocks in the Earth's mantle or rocksalt in the salt diapirs, is nearly incompressible, which is well-known to require a stable discretization to avoid numerical artifacts, that include volumetric and shear locking, pressure oscillations, and hourglass modes (e.g., Zienkiewicz and Taylor, 2000). Moreover, the same numerical problems occur when plastic deformation dominates irrespective whether the flow is incompressible or dilatant/contractant, since the volumetric plastic deformation can be effectively viewed as a form of kinematical constraint, which is similar to incompressibility (e.g., de Borst and Groen, 1995).

Within the finite element framework, a typical approach to deal with incompressibility is to use a so-called mixed formulation which combines velocity (displacement increments) and pressure as primary variables. Yet not all interpolation types can be adopted since they need to satisfy a set of stability criteria commonly known as the LBB conditions (e.g., Gatica, 2014). Examples of stable finite elements types include the Taylor-Hood approximation on quadrilateral and hexagonal shape e.g. Q₂×Q₁, which uses continuous quadratic and linear polynomials for the velocity and pressure, respectively. The discontinuous pressure version is typically represented by $Q_{\mathrm{2}}\times P_{-\mathrm{1}}$ elements (see e.g., Thieulot and Bangerth, 2022, for an overview of most commonly used combinations). The latter element is very widely used in geodynamic simulations due to its reliability in dealing with sharp jumps in viscosities (Kronbichler et al., 2012; May et al., 2014; Deubelbeiss and Kaus, 2008).

The simplex element shape can also be used to construct stable discretizations. A particular candidate is the conforming triangular Crouzeix-Raviart element ($P_{\mathrm{2}}^{+}\times P_{-\mathrm{1}}$), which uses quadratic interpolation enhanced with a bubble function for the velocities (displacement increments), and a linear discontinuous interpolation for pressure (Crouzeix and Raviart, 1973). The triangular shape of the element facilitates discretizing complex geometrical domains, whereas the discontinuous pressure interpolation enables enforcing mass conservation at the element level and is in general advantageous for difficult problems with abrupt material contrasts (Pelletier et al., 1989). The element has been shown to behave robustly in practice (Dabrowski et al., 2008), and can be generalized to 3D (Crouzeix and Raviart, 1973).

As an alternative to the stable finite elements, it is also possible to apply a staggered grid finite difference discretization (Harlow and Welch, 1965). This method is computationally cheap and also demonstrates excellent and reliable performance for the relevant problems (Gerya and Yuen, 2007; Tackley, 2008; Kaus et al., 2016; Räss et al., 2017; Deubelbeiss and Kaus, 2008). The staggered finite difference formulation was proven to be a stable formulation for the incompressible (Shin and Strikwerda, 1997), and compressible Stokes equations (Eymard et al., 2010), and was numerically demonstrated to be LBB-stable (Gerya et al., 2013).

Despite that a number of stabilized equal-order interpolation finite elements have been proposed (e.g. Dohrmann and Bochev, 2004; Cioncolini and Boffi, 2019), their performance does not practically bring a satisfactory level of robustness for the cases of large parameter variation (Thieulot and Bangerth, 2022). Selecting a stable discretization is therefore an important milestone for successful numerical implementation of the nonlinear rheological models that involve nearly-incompressible material behavior.

Implementing plastic rheologies in the context of a two-field formulation is not straightforward. The problem is caused by the presence of stress-like variables (pressure) among the primary unknowns, whereas stress integration algorithms are typically formulated in a strain-driven fashion (de Borst et al., 2012). Pressure-independent models, such as J2-plasticity (Pastor et al., 1997), or incompressible Drucker-Prager plasticity (Gerya and Yuen, 2007; May et al., 2014; Kaus et al., 2016; Glerum et al., 2018) poses no additional difficulty because pressure remains unaltered during the local stress update. For nonzero dilatation cases, as well as for plastic compaction or tensile failure modes, this is no longer the case, and therefore at least two potential solutions may be elaborated. The first approach computes deviatoric stresses from the kinematical variables (velocities or displacement increments) using the standard strain-driven algorithm, wheres the globally discretized pressure is treated as an actual spherical part of the Cauchy stress tensor (Commend, 2001; Commend et al., 2004). This type of algorithms can be referred to as true pressure formulation. The second approach treats the globally discretized pressure as a trial visco-elastic pressure (Duretz et al., 2021), which can be termed as trial pressure formulation. Additionally, it should be noted that pressure selection problems can be completely avoided by adopting a mixed formulation that involves strains as global variables (e.g., Benedetti et al., 2014). This approach is computationally much more expensive and therefore beyond the scope of this paper.

It is often quite difficult to describe different plastic failure regimes, such as plastic compaction, or tensile failure (mode I) and shear failure (mode II) with a single yield surface. A standard approach is therefore to tackle this problem with a multi-surface plasticity model. A two- or three-invariant shear failure envelope is typically combined with compression and tensile caps (e.g., Sandler and Rubin, 1979). An example of a plasticity model that handles a combination of tensile and shear failure of concrete in plane stress formulation can be found in Feenstra and de Borst (1996). In general, a non-smooth intersection of the yield surface segments results in the complex algorithmic treatment of the corner regions in the stress space (e.g., Simo et al., 1988). To avoid this difficulty, alternative smooth multi-surface models have been proposed where the focus was mostly placed on enforcing continuity and differentiability of the yield surface (Swan and Seo, 2000; Dolarevic and Ibrahimbegovic, 2007; Motamedi and Foster, 2015). This, however, does not automatically guarantee continuity in the entire stress space. Plotting isocontours of the yield surface can reveal typical problems that include spurious elastic domains, singularity points, erratic gradients, discontinuities, or non-convex regions (e.g., Stupkiewicz et al., 2014; Golchin et al., 2021). A squared version of the yield functions does not even ensure dimensional consistency between the segments (Swan and Seo, 2000), or may lead to loss of convexity in the stress space (Stupkiewicz et al., 2014). Complex stress integration procedures are therefore required to determine the set of active yield surface segments (Swan and Seo, 2000; Dolarevic and Ibrahimbegovic, 2007). All these problems lead to severe algorithmic difficulties and convergence problems of local stress update iterations.

To model rock failure, the tensile regime is more important than plastic compaction as it is relevant for modeling dyke propagation (e.g., Rivalta et al., 2015; Li et al., 2023). It is therefore potentially advantageous to describe both failure modes with a single smooth hyperbolic approximation of the Mohr-Coulomb yield surface (Abbo and Sloan, 1995). This approach, however, completely fails when dilatation angles go to zero, since the plastic flow potential becomes pressure-insensitive in the tensile mode. Carol et al. (1997) proposed an ad-hoc approach to restore the proper tensile behavior by introducing a pressure-dependent dilation coefficient, however, at the cost of destroying the convexity of the flow potential, which can be easily revealed by plotting its isocontours. Li et al. (2023) adopted the same approximation, but in their formulation, the amount of dilatation in flow potential is dependent on porosity rather than on pressure, which lacks physical explanation. Other existing implementations of tensile failure of rocks in the geodynamic community (e.g., Rozhko et al., 2007; Keller et al., 2013; Kiss et al., 2023) employ non-smooth yield surfaces, and are therefore not guaranteed to produce well-defined stress integration algorithms at least in the context of an implicit time integration scheme.

Classical rate-independent plasticity models that incorporate strain softening, which is typically required by constitutive models of rocks and soils (e.g., Read and Hegemier, 1984), are mathematically ill-posed. This results in the lack of an intrinsic length-scale, severe mesh-dependence and convergence issues (de Borst and Mühlhaus, 1992; Spiegelman et al., 2016). There are several approaches to remediate this problem (e.g. Duretz et al., 2023), with the easiest one being a viscoplastic rate-dependent regularization (de Borst and Duretz, 2020). The latter assumes different ways to formulate the constitutive equation for the viscoplastic strain rate. Namely, Perzyna model (Perzyna, 1966) assumes an explicit equation, Duvaut-Lions model (Duvaut and Lions, 1972) is built upon a relatively simple generalization of a rate-independent plasticity, whereas consistency model (Wang et al., 1997; Heeres et al., 2002) introduces the concept of a rate-dependent yield function. Irrespective of the employed viscoplastic regularization, either Perzyna (e.g., Jacquey and Cacace, 2020a, b), Duvaut-Lions (e.g., Kiss et al., 2023), or the consistency model (e.g., Duretz et al., 2019; Li et al., 2023), a successful resolution of the mesh issue simultaneously with the improvement of the global convergence or stability of the time integration scheme can be achieved.

To avoid potential misconceptions, it should be noted that despite satisfactory practical performance of viscoplastic regularizations, a rigorous theoretical framework proving its effectiveness in quasi-static cases, such as considered in this work, does not currently exist (e.g., de Borst and Duretz, 2020). Moreover, even for dynamic cases, questions have recently been raised on the regularization potential of viscoplasticity, with some suggesting that it only works conditionally (e.g., Jacquey et al., 2021), whereas others suggest that it is unable to regularize a 1D simple shear setup (e.g., Stathas and Stefanou, 2022). Despite the lack of a solid theoretical background, we still advocate the use of viscoplastic regularization motivated by its successful application to multi-dimensional quasi-static cases, which can be traced back to the earlier works of Zienkiewicz and Cormeau (1972) and Zienkiewicz and Cormeau (1974) and which is fully supported by more recent publications we cite above, as well as the results shown here. It should also be mentioned that viscoplasticity extends beyond the regularization framework and may represent a true physical deformation mechanism that requires laboratory experiments to calibrate its material parameters.

Building up on previous work, we present a visco-elasto-viscoplastic rheological model that combines linear elastic response with diffusion and dislocation viscous creep mechanisms and a relatively simple Perzyna-type viscoplastic model. The yield surface is composed of a linear Drucker-Prager shear failure envelop and a circular tensile cap in a way that ensures dimensional consistency, convexity, continuity and differentiability throughout the entire stress space. The model allows an arbitrary amount of dilatation in the shear failure regime, without compromising the description of tensile failure. The viscoplastic regularization enables the incorporation of strain softening, such that spurious mesh dependence is avoided and global convergence is ensured. We provide algorithmic details necessary to compute the stresses at the integration points and to solve the resulting global system of nonlinear equilibrium equations with a Newton-Raphson method, including the tangent matrix derivation. The algorithms are implemented in a 2D finite element code that employs stable conforming Crouzeix-Raviart elements, along with incremental displacements and incremental trial pressure as the primary unknowns. A number of examples are given to demonstrate the numerical robustness of the code. The plasticity model presented in this paper can be readily implemented in any finite element or finite difference code that uses a stable mixed two-field formulation.

2 Physical model

2.1 Volumetric-deviatoric decomposition

Throughout this paper, we assume a standard decomposition of an arbitrary tensor (e.g. a_ij) into its volumetric (spherical) and deviatoric projections, which are given by, respectively, in the index notation (Einstein summation convention is implied):

$$\begin{array}{}\text{(1)}& \mathrm{tr}\left(a_{ij}\right)=a_{kk},\mathrm{dev}\left(a_{ij}\right)=a_{ij}-{\displaystyle \frac{\mathrm{1}}{\mathrm{3}}}{a}_{kk}{\mathit{\delta }}_{ij},\end{array}$$

where δ_ij is the Kronecker delta or second order unit tensor: δ_ij=1, for i=j, δ_ij=0, for i≠j. For an arbitrary deviatoric tensor we additionally introduce an effective scalar measure equal to the square root of its second invariant, defined as:

$$\begin{array}{}\text{(2)}& a_{\mathrm{II}}={\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{a}_{ij}{a}_{ij}\right)}^{\frac{\mathrm{1}}{\mathrm{2}}}.\end{array}$$

For the Cauchy stress tensor (σ_ij) the volumetric-deviatoric decomposition into stress deviator (τ_ij) and mean stress or pressure (p), which is assumed to be positive in compression, is defined as:

$$\begin{array}{}\text{(3)}& {\mathit{\sigma }}_{ij}={\mathit{\tau }}_{ij}-p{\mathit{\delta }}_{ij},{\mathit{\tau }}_{ij}=\mathrm{dev}\left({\mathit{\sigma }}_{ij}\right),p=-{\displaystyle \frac{\mathrm{1}}{\mathrm{3}}}\mathrm{tr}\left({\mathit{\sigma }}_{ij}\right).\end{array}$$

Equivalently, the strain rate tensor (${\dot{\mathit{ϵ}}}_{ij}$) can be decomposed into a strain rate deviator (${\dot{\mathit{\epsilon }}}_{ij}$) and a volumetric strain rate ($\dot{\mathit{\theta }}$) by:

$$\begin{array}{}\text{(4)}& \begin{array}{rl}& {\dot{\mathit{ϵ}}}_{ij}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left({\displaystyle \frac{\partial v_i}{\partial x_j}}+{\displaystyle \frac{\partial v_j}{\partial x_i}}\right)={\dot{\mathit{\epsilon }}}_{ij}+{\displaystyle \frac{\mathrm{1}}{\mathrm{3}}}\dot{\mathit{\theta }}{\mathit{\delta }}_{ij},\\ & {\dot{\mathit{\epsilon }}}_{ij}=\mathrm{dev}\left({\dot{\mathit{ϵ}}}_{ij}\right),\\ & \dot{\mathit{\theta }}=\mathrm{tr}\left({\dot{\mathit{ϵ}}}_{ij}\right),\end{array}\end{array}$$

where $x_i,i=\mathrm{1},\mathrm{2},\mathrm{3}$ denote Cartesian coordinates, and v_i are the components of the spatial velocity vector.

2.2 Constitutive and conservation equations

The constitutive equations for the deviatoric stress assumes additive decomposition of the deviatoric strain rate into elastic, diffusion creep, dislocation creep and viscoplastic components, respectively, as follows:

$$\begin{array}{}\text{(5)}& \begin{array}{rl}{\dot{\mathit{\epsilon }}}_{ij}=& {\dot{\mathit{\epsilon }}}_{ij}^{\mathrm{el}}+{\dot{\mathit{\epsilon }}}_{ij}^{\mathrm{dif}}+{\dot{\mathit{\epsilon }}}_{ij}^{\mathrm{dis}}+{\dot{\mathit{\epsilon }}}_{ij}^{\mathrm{vp}}={\displaystyle \frac{{\stackrel{\diamond }{\mathit{\tau }}}_{ij}}{\mathrm{2}G}}+A_{\mathrm{D}}{\mathit{\tau }}_{ij}\\ & +A_{\mathrm{N}}{\left({\mathit{\tau }}_{\mathrm{II}}\right)}^{n-\mathrm{1}}{\mathit{\tau }}_{ij}+\dot{\mathit{\lambda }}\mathrm{dev}\left({\displaystyle \frac{\partial Q}{\partial {\mathit{\sigma }}_{ij}}}\right).\end{array}\end{array}$$

Here G denotes the elastic shear modulus, A_D and A_N are the diffusion and dislocation creep prefactors, respectively, n is the power-law exponent, Q is the plastic flow potential function, and $\dot{\mathit{\lambda }}$ is the magnitude of the viscoplastic strain rate (viscoplastic multiplier), which can be explicitly formulated for the Perzyna model (see Sect. 2.3 and 2.4). The Jaumann objective stress rate is defined as:

$$\begin{array}{}\text{(6)}& {\stackrel{\diamond }{\mathit{\tau }}}_{ij}={\displaystyle \frac{\partial {\mathit{\tau }}_{ij}}{\partial t}}+{\mathit{\tau }}_{ik}{\mathit{\omega }}_{kj}-{\mathit{\omega }}_{ik}{\mathit{\tau }}_{kj},\end{array}$$

where ${\dot{\mathit{\omega }}}_{ij}$ is the spin tensor, given by:

$$\begin{array}{}\text{(7)}& {\dot{\mathit{\omega }}}_{ij}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left({\displaystyle \frac{\partial v_i}{\partial x_j}}-{\displaystyle \frac{\partial v_j}{\partial x_i}}\right).\end{array}$$

The temperature-dependent creep prefactors are:

$$\begin{array}{}\text{(8)}& A_{\mathrm{D}}=B_{\mathrm{D}}\exp \left[-{\displaystyle \frac{E_{\mathrm{L}}}{RT}}\right],A_{\mathrm{N}}=B_{\mathrm{N}}\exp \left[-{\displaystyle \frac{E_{\mathrm{N}}}{RT}}\right],\end{array}$$

where B_D, B_N and E_D, E_N denote the corresponding constants and activation enthalpies of the diffusion and dislocation creep, respectively, T is the temperature, and R the gas constant.

The volumetric constitutive equation, which can also be termed as continuity equation, is defined similarly:

$$\begin{array}{}\text{(9)}& \dot{\mathit{\theta }}={\dot{\mathit{\theta }}}^{\mathrm{el}}+{\dot{\mathit{\theta }}}^{\mathrm{vp}}=-{\displaystyle \frac{\mathrm{1}}{K}}{\displaystyle \frac{\mathrm{D}p}{\mathrm{D}t}}+\dot{\mathit{\lambda }}\mathrm{tr}\left({\displaystyle \frac{\partial Q}{\partial {\mathit{\sigma }}_{ij}}}\right).\end{array}$$

Here the total volumetric strain rate is additively decomposed into an elastic and viscoplastic component, K is the elastic bulk modulus, and $\mathrm{D}/\mathrm{D}t$ is the material time derivative. A simple uni-axial idealization of both deviatoric and volumetric constitutive equations is illustrated in Fig. 1.

Figure 1 Schematic uni-axial representation of the rheological model. Elastic, viscous and plastic deformations are idealized by the spring, dashpot, and sliding frictional elements, receptively. (a) Deviatoric constitutive equation. (b) Volumetric constitutive equation.

In case the volumetric response becomes much stiffer compared to the deviatoric response (quasi-incompressible limit), it becomes mandatory to revert to a mixed formulation to compute accurate pressure fields. The continuity equation is discretized at the global level, and pressure is treated as a global independent variable. For incompressible Stokes flow, the continuity equation simply reduces to a constraint of the form $\dot{\mathit{\theta }}=\mathrm{0}$, and pressure becomes a corresponding Lagrange multiplier variable. Finally, the computed deviatoric stresses and pressure should satisfy the global equilibrium equation also known as Cauchy momentum equation, given by:

$$\begin{array}{}\text{(10)}& {\displaystyle \frac{\partial {\mathit{\tau }}_{ij}}{\partial x_j}}-{\displaystyle \frac{\partial p}{\partial x_i}}+\mathit{\rho }{g}_i=\mathrm{0},\end{array}$$

where ρ is the material density, and g_i are the components of the gravity acceleration vector. Here, we only consider quasi-static processes and hence explicitly ignore the inertial terms in the momentum equation.

2.3 Perzyna viscoplasticity

To deal with mesh dependency issues and convergence problems, which are typical for rate-independent plasticity models (e.g., Spiegelman et al., 2016; Duretz et al., 2019) we adopt a viscoplastic formulation following Perzyna (1966). In contrast with the Duvaut-Lions or consistency models, it has an explicit equation that defines the components of the viscoplastic strain rate. Here we use the simplest form of the Perzyna model without the power-law function, as it is defined in e.g. Jacquey and Cacace (2020a, b). The magnitude of the viscoplastic multiplier can be written as:

$$\begin{array}{}\text{(11)}& \dot{\mathit{\lambda }}={\displaystyle \frac{\langle F\rangle }{{\mathit{\eta }}^{\mathrm{vp}}}},\end{array}$$

where η^vp is the viscoplastic regularization viscosity, F is the rate-independent plastic yield function surrounded by Macaulay brackets, which has the following meaning:

$$\begin{array}{}\text{(12)}& \langle F\rangle =\left\{\begin{array}{ll}F,& F\ge \mathrm{0}\\ \mathrm{0},& \mathrm{otherwise}\end{array}\right..\end{array}$$

The rate-independent plasticity is naturally recovered with this formulation for η^vp→0. Likewise, the viscoplastic deformation can be effectively switched off and the visco-elastic solution can be enforced with the limit η^vp→∞.

Irrespective of the particular viscoplastic model, an addition of a regularization term implicitly introduces a length scale in strain localization problems (e.g. Duretz et al., 2023). It is important to note that regularization viscosity should not be treated as a part of the rheological model, but rather as a numerical parameter, which requires careful selection depending on the process being modeled (Duretz et al., 2019). There is also a trade-off between the sharpness of localization and number of nonlinear iterations required to achieve global convergence. In other words, increasing the amount of regularization improves the convergence but simultaneously smears out the localization zone. In the context of the Perzyna model, another issue requires additional attention. Direct inspection of Eq. (11) reveals that viscoplastic strain initiates only if stress violates the yield function constraint, i.e. a certain amount of overstress occurs. This feature imposes an explicit restriction on the yield (overstress) function to be continuous and convex in the entire stress space (Simo, 1989).

2.4 Smooth yield surface

Here we propose a simple two-surface plasticity model that combines a linear Drucker-Prager shear failure envelop with a circular tensile cap function in a manner that ensures its applicability in the context of Perzyna viscoplasticity (Fig. 2b). A similar model was proposed by Swan and Seo (2000), but continuity was only enforced at the yield surface, and dimensional consistency was violated due to usage of the squared versions of the yield functions. Our model requires four independent material parameters, which include the friction angle (φ), the dilation angle (ψ), the Mohr-Coulomb cohesion (c_MC), and the tensile strength (p_T). The other model parameters can be computed from the primary ones via geometrical transformations, which are omitted here for clarity. The overall parameter layout is schematically illustrated in Fig. 2d.

Figure 2 Meridional plot of the smooth yield function and flow potential. (a) Yield function map without scaling (a=1 in Eq. 18). (b) Yield function map with scaling (a≠1 in Eq. 18). (c) Flow potential map above the yield surface (F>0). (d) Schematic illustration of the yield surface and flow potential parameters. Refer to the code availability section to access the script that generates the yield surface and flow potential maps.

The Drucker-Prager cohesion (c), friction coefficient (k) and dilatation coefficient (k_q) are computed as follows:

$$\begin{array}{}\text{(13)}& k=\sin \mathit{\phi },k_q=\sin \mathit{\psi },c=c_{\mathrm{MC}}\cos \mathit{\phi }.\end{array}$$

This approach is rather simple, widely used in geodynamics (Keller et al., 2013; Kaus et al., 2016; Kaus, 2010; Duretz et al., 2021), and ensures that frictional strength is not systematically overestimated. For the details and alternative methods we refer to e.g. Jiang and Xie (2011). To ensure that the composite yield surface produces a continuous and differential map two issues need to be addressed: (i) a delimiter point between the segments must be selected such that they intersect in a smooth manner (Fig. 2a), and (ii) the tensile cap function must be scaled to eliminate the discontinuity outside the yield surface (Fig. 2b). We therefore introduce the following scaling coefficients, one for the yield surface (a), and one for the flow potential (b):

$$\begin{array}{}\text{(14)}& a=\sqrt{\mathrm{1}+k^{\mathrm{2}}},b=\sqrt{\mathrm{1}+k_q^{\mathrm{2}}}.\end{array}$$

The center coordinate and the radius of the cap function, respectively, can be computed by enforcing the smooth intersection between both segments as:

$$\begin{array}{}\text{(15)}& p_{\mathrm{y}}=\left(p_{\mathrm{T}}+{\displaystyle \frac{c}{a}}\right){\left(\mathrm{1}-{\displaystyle \frac{k}{a}}\right)}^{-\mathrm{1}},R_{\mathrm{y}}=p_{\mathrm{y}}-p_{\mathrm{T}}.\end{array}$$

The coordinates of the delimiter point are therefore given by:

$$\begin{array}{}\text{(16)}& p_{\mathrm{d}}=p_{\mathrm{y}}-R_{\mathrm{y}}{\displaystyle \frac{k}{a}},{\mathit{\tau }}_{\mathrm{d}}=k{p}_{\mathrm{d}}+c.\end{array}$$

Finally, the center coordinate of the cap flow potential can be assigned such that the transition between the tensile and shear domains in the stress space passes through the delimiter point (see Fig. 2c), which yields:

$$\begin{array}{}\text{(17)}& p_q=p_{\mathrm{d}}+k_q{\mathit{\tau }}_{\mathrm{d}}.\end{array}$$

Having defined all necessary model parameters we are ready to write down the equations for the yield function:

$$\begin{array}{}\text{(18)}& F=\left\{\begin{array}{ll}{\mathit{\tau }}_{\mathrm{II}}-kp-c,& {\mathit{\tau }}_{\mathrm{II}}\left(p_{\mathrm{y}}-p_{\mathrm{d}}\right)\ge {\mathit{\tau }}_{\mathrm{d}}\left(p_{\mathrm{y}}-p\right)\\ a\left({\widehat{R}}_y-R_{\mathrm{y}}\right),& \mathrm{otherwise}\end{array}\right.,\end{array}$$

and the corresponding flow potential:

$$\begin{array}{}\text{(19)}& Q=\left\{\begin{array}{ll}{\mathit{\tau }}_{\mathrm{II}}-k_{q}p-\mathrm{const},& {\mathit{\tau }}_{\mathrm{II}}\left(p_q-p_{\mathrm{d}}\right)\ge {\mathit{\tau }}_{\mathrm{d}}\left(p_q-p\right)\\ b\left({\widehat{R}}_q-\mathrm{const}\right),& \mathrm{otherwise}\end{array}\right..\end{array}$$

Here the radii of the points in the stress space for both the tensile yield function and flow potential, respectively, are given by:

$$\begin{array}{}\text{(20)}& {\widehat{R}}_y=\sqrt{{\mathit{\tau }}_{\mathrm{II}}^{\mathrm{2}}+{\left(p-p_{\mathrm{y}}\right)}^{\mathrm{2}}},{\widehat{R}}_q=\sqrt{{\mathit{\tau }}_{\mathrm{II}}^{\mathrm{2}}+{\left(p-p_q\right)}^{\mathrm{2}}}.\end{array}$$

It should be pointed out that our model gives a simple analytical description of the stress space partitioning into shear and tensile domains, and therefore no complex algorithm is required to detect the active yield surface. In any point above the yield surface the gradient of the flow potential is uniquely defined, and hence the return mapping direction is completely constrained. This property greatly simplifies the implementation of the stress update algorithms. Since an arbitrary amount of dilatation is allowed in our model, the resulting flow potential function is non-associated, as is illustrated in Fig. 2c. The dilatation angle for the Drucker-Prager part of the yield function does not affect the description of the tensile failure, which implies that zero dilatation angle cases are explicitly supported.

We complete the formulation of the plasticity model by directly differentiating Eq. (19) to obtain the required expressions for the flow potential gradients:

$$\begin{array}{}\text{(21)}& {\displaystyle \frac{\partial Q}{\partial {\mathit{\sigma }}_{ij}}}=B_{\mathit{\tau }}{\mathit{\tau }}_{ij}+B_{\mathrm{p}}{\mathit{\delta }}_{ij},\end{array}$$

where the prefactors are given by:

$$\begin{array}{}\text{(22)}& \begin{array}{rl}& \left(B_{\mathit{\tau }},B_{\mathrm{p}}\right)=\\ & \left\{\begin{array}{ll}\left(\frac{\mathrm{1}}{\mathrm{2}{\mathit{\tau }}_{\mathrm{II}}},\frac{k_q}{\mathrm{3}}\right),& {\mathit{\tau }}_{\mathrm{II}}\left(p_q-p_{\mathrm{d}}\right)\ge {\mathit{\tau }}_{\mathrm{d}}\left(p_q-p\right)\\ \left(\frac{b}{\mathrm{2}{\widehat{R}}_q},-\frac{b(p-p_q)}{\mathrm{3}{\widehat{R}}_q}\right),& \mathrm{otherwise}\end{array}\right..\end{array}\end{array}$$

2.5 Strain softening

Plastic deformation in the rocks and soils is generally characterized by the occurrence of the strain softening (see the extensive review by Read and Hegemier, 1984, for the physical details and explanations of this phenomenon) . Here, we parametrize strain softening for the Mohr-Coulomb friction angle and cohesion as a linear function of the accumulated deviatoric viscoplastic strain, defined as:

$$\begin{array}{}\text{(23)}& \mathit{\kappa }=\underset{t}{\int }{\dot{\mathit{\epsilon }}}_{\mathrm{II}}^{\mathrm{vp}}\mathrm{d}t.\end{array}$$

The expressions for the friction angle and cohesion can be formulated as follows:

$$\begin{array}{}\text{(24)}& \begin{array}{rl}& \mathit{\phi }=\max \left({\mathit{\phi }}^{\mathrm{init}}+H_{\mathit{\phi }}\mathit{\kappa },{\mathit{\phi }}^{\min }\right),\\ & c_{\mathrm{MC}}=\max \left(c_{\mathrm{MC}}^{\mathrm{init}}+H_{\mathrm{c}}\mathit{\kappa },c_{\mathrm{MC}}^{\min }\right),\end{array}\end{array}$$

where the superscripts ^init and ^min denote the initial and saturated values for the corresponding material parameters, H_φ and H_c are the hardening modulii for the friction angle and cohesion, respectively. Note that the hardening modulii are negative for softening cases.

To identify zones of tensile failure and dilatant shear plasticity the accumulated volumetric viscoplastic strain can be defined in a similar way:

$$\begin{array}{}\text{(25)}& \mathit{\chi }=\underset{t}{\int }{\dot{\mathit{\theta }}}^{\mathrm{vp}}\mathrm{d}t.\end{array}$$

3 Numerical formulation

3.1 Time discretization

To integrate the coupled nonlinear constitutive equations in time we apply backward Euler implicit time discretization which is unconditionally stable and first order accurate. All rate quantities are presented by their unknown values at the end of the time step and assumed to be constant during the time step. In that case the computation of the incremental quantities can be done trivially, e.g. for the deviatoric strain rate tensor we can write:

$$\begin{array}{}\text{(26)}& \mathrm{\Delta }{\mathit{\epsilon }}_{ij}={\dot{\mathit{\epsilon }}}_{ij}\mathrm{\Delta }t,\end{array}$$

where $\mathrm{\Delta }t=t_{n+\mathrm{1}}-t_n$ is the time step, t_n is the time in the beginning, and t_n+1 is the time in the end of the nth time step. Similarly, the rate quantities can be estimated from their corresponding increments as:

$$\begin{array}{}\text{(27)}& {\dot{\mathit{\epsilon }}}_{ij}={\displaystyle \frac{\mathrm{\Delta }{\mathit{\epsilon }}_{ij}}{\mathrm{\Delta }t}}.\end{array}$$

To simplify the notation we omit the indices at the end of the time step for all quantities in the remainder of this paper, e.g. for the pressure we write p instead of p_n+1. For the converged values from the previous time step we always include the time step index, e.g. for the history pressure we write p_n.

3.2 True pressure scheme

In the context of a mixed formulation, the globally discretized velocity (displacement increment) and pressure are typically treated as primary unknowns. Since dilatant plasticity, in general, involves pressure modifications during local stress updates, we must decide how to interpret the global pressure variable (which is obtained from solving the global system of nonlinear equations). One of the approaches is based on treating the global pressure as a true spherical part of the Cauchy stress tensor (Commend et al., 2004). Hence, the deviatoric stress, local pressure, and viscoplastic volumetric strain rate must be computed using a standard strain-driven approach, solely as the functions of the strain rate, and hence the velocity:

$$\begin{array}{}\text{(28)}& {\mathit{\tau }}_{ij}={\mathit{\tau }}_{ij}\left({\dot{\mathit{ϵ}}}_{ij}\right),p=p\left({\dot{\mathit{ϵ}}}_{ij}\right),{\dot{\mathit{\theta }}}^{\mathrm{vp}}={\dot{\mathit{\theta }}}^{\mathrm{vp}}\left({\dot{\mathit{ϵ}}}_{ij}\right).\end{array}$$

At this stage, the globally discretized pressure, denoted as p^*, may differ from the locally computed pressure in the integration points, i.e.: $p^{*}\ne p$. Since the global pressure is treated as a true pressure variable, the local pressure should be discarded, and only the deviatoric stress and viscoplastic volumetric strain rate should be used for the evaluation of the time-discrete momentum and continuity residual equations, respectively:

$$\begin{array}{}\text{(29)}& \begin{array}{rl}& r_i^{\mathrm{m}}={\displaystyle \frac{\partial {\mathit{\tau }}_{ij}}{\partial x_j}}-{\displaystyle \frac{\partial p^{*}}{\partial x_i}}+\mathit{\rho }{g}_i,\\ & r^{\mathrm{c}}=\dot{\mathit{\theta }}+{\displaystyle \frac{\mathrm{1}}{K}}{\displaystyle \frac{p^{*}-p_n}{\mathrm{\Delta }t}}-{\dot{\mathit{\theta }}}^{\mathrm{vp}}.\end{array}\end{array}$$

The addition of the viscoplastic strain rate to the continuity residual ensures that the global pressure receives a feedback from the volumetric viscoplasticity. Upon convergence of the global nonlinear iterations, both pressure variables should become approximately equal, i.e.: $p^{*}\approx p$.

3.3 Trial pressure scheme

The globally discretized pressure variable (p^*) can alternatively be interpreted as a trial visco-elastic pressure (Duretz et al., 2021). In this case, the deviatoric stress, local pressure, and viscoplastic volumetric strain rate become functions of both strain rate and global pressure:

$$\begin{array}{}\text{(30)}& \begin{array}{rl}& {\mathit{\tau }}_{ij}={\mathit{\tau }}_{ij}\left({\dot{\mathit{ϵ}}}_{ij},p^{*}\right),p=p\left({\dot{\mathit{ϵ}}}_{ij},p^{*}\right),\\ & {\dot{\mathit{\theta }}}^{\mathrm{vp}}={\dot{\mathit{\theta }}}^{\mathrm{vp}}\left({\dot{\mathit{ϵ}}}_{ij},p^{*}\right).\end{array}\end{array}$$

Specifically, the pressure update assumes the following discrete form in this formulation (Duretz et al., 2021):

$$\begin{array}{}\text{(31)}& p=p^{*}+K{\dot{\mathit{\theta }}}^{\mathrm{vp}}\mathrm{\Delta }t.\end{array}$$

We first formulate the time-discrete residual equations involving both the updated deviatoric stress and pressure:

$$\begin{array}{}\text{(32)}& {\displaystyle }{r}_i^{\mathrm{m}}={\displaystyle \frac{\partial {\mathit{\tau }}_{ij}}{\partial x_j}}-{\displaystyle \frac{\partial p}{\partial x_i}}+\mathit{\rho }{g}_i,\text{(33)}& {\displaystyle }{r}^{\mathrm{c}}=\dot{\mathit{\theta }}+{\displaystyle \frac{\mathrm{1}}{K}}{\displaystyle \frac{p-p_n}{\mathrm{\Delta }t}}-{\dot{\mathit{\theta }}}^{\mathrm{vp}}.\end{array}$$

Substituting Eq. (31) into the continuity residual equation (Eq. 33) yields:

$$\begin{array}{}\text{(34)}& r^{\mathrm{c}}=\dot{\mathit{\theta }}+{\displaystyle \frac{\mathrm{1}}{K}}{\displaystyle \frac{p^{*}-p_n}{\mathrm{\Delta }t}}.\end{array}$$

The viscoplastic volumetric strain rate simply cancels from the global continuity residual. The reason for this is that the local pressure update (Eq. 31) is derived by enforcing the continuity equation in its complete form (Eq. 33). This implies that the continuity residual (Eq. 33) is always balanced at the level of the local stress update. The local pressure (p) represents the true spherical part of the Cauchy stress tensor upon achieving global convergence. The difference between the global and local pressure does not vanish, i.e.: $p^{*}\ne p$. The global pressure (p^*) converges to the trial pressure value. The global solver should therefore directly apply spatial discretization to the Eqs. (32) and (34). Using the Eq. (33) in the global solver is not justified, since it is indistinguishable from the Eq. (34).

3.4 Pressure scheme comparison

Both pressure schemes should theoretically deliver the same results. However their convergence properties are not guaranteed to be the same. First, it should be noted that for non-dilatant plasticity cases, the difference between the pressure variables does not occur in either of the approaches, i.e. $p^{*}=p$. It might also be taken for granted that both schemes should deliver the same linearization in this case. This assumption does not hold when the yield function still depends on pressure, i.e. F=F(p), such as for Drucker-Prager plasticity with a zero dilatation coefficient. The true pressure scheme would still imply that pressure used to evaluate the yield function is estimated from the strain rates, i.e.: $p=p\left({\dot{\mathit{ϵ}}}_{ij}\right)$. The trial pressure scheme would instead directly use the global pressure (p^*). Even when both pressure values should theoretically be the same, at least for stable discretizations, this still formally leads to different functional dependencies, different linearizations, and different convergence properties. Only for the truly pressure-independent plasticity models, like the Von Mises model, both schemes are guaranteed to be exactly same (which is irrelevant in the context of rock failure modeling). It should be also noted that both pressure schemes always deliver non-symmetric global tangent matrices for the pressure-dependent plasticity models. This statement holds even for completely associated models, when the friction angle is equal to the dilatation angle. In general, all dilatant plasticity models require the elastic bulk modulus (K) to be finite, and both pressure schemes will fail to handle dilatant plasticity cases for elastically incompressible materials (K→∞).

A detailed analysis of the convergence properties of both pressure schemes goes beyond the scope of this paper. Here, we simply report our practical observations and conclusions. We have thoroughly implemented, linearized and tested both pressure schemes for the relevant dilatant and non-dilatant plasticity cases (see Berlie, 2023, for further details). The results demonstrate that the trial pressure scheme performs more robustly compared to the true pressure formulation. In the following, we therefore limit our discussion to the trial pressure scheme.

3.5 Primary variable selection

Traditionally, velocity is selected as the primary kinematical variable to solve boundary value problems that involve quasi-static deformation of viscous materials. However, this is a suboptimal choice in the presence of highly nonlinear rheologies and Dirichlet boundary conditions. The problem is caused by the lack of time step control over the magnitude of the kinematical loading. With a velocity formulation, the entire boundary velocity is applied instantly, which may result in severe convergence problems. There is essentially no measure that can be applied to mitigate this problem, since boundary velocities are independent of the time step. We therefore choose to use displacement increments (Δu_i) rather than velocity as the primary kinematical variables in this work. This approach has the advantage that kinematical loading can be directly controlled by the time step, since the velocity and displacement increments are related in the following manner:

$$\begin{array}{}\text{(35)}& \mathrm{\Delta }{u}_i=v_i\mathrm{\Delta }t.\end{array}$$

This property opens the possibility to design an adaptive time stepping algorithm. As soon as the selected time step is too large to achieve convergence within a reasonable number of iterations, the equilibrium iteration can be simply restarted with a smaller time step magnitude, which implies a smaller displacement step and therefore a smaller stress increment.

We emphasize that using displacement increments instead of velocities imposes no restrictions on the type of rheology that can be handled (viscous, elastic, viscoplastic, or combinations thereof). Both velocities and strain rates are perfectly defined in this formulation (see Eq. 27). Hence the stress integration algorithm can be formulated using strain rates. For instance the algorithm presented in Sect. 3.6 can be readily used in the framework of a finite element or a finite difference code that uses a standard approach with the velocity as a primary variable.

As the second global primary variable in the mixed formulation, we select the trial pressure increment (Δp^*). The total trial pressure is thus defined as:

$$\begin{array}{}\text{(36)}& p^{*}=p_n+\mathrm{\Delta }{p}^{*}.\end{array}$$

3.6 Local stress update algorithm

The integration of constitutive equations that include a viscoplastic component is usually done with a two-stage predictor-corrector procedure (see Appendix A for derivations). During the first stage, the magnitude of the viscoplastic multiplier ($\dot{\mathit{\lambda }}$) is assumed to be zero and the entire strain rate is visco-elastic (predictor stage). The magnitude of the trial visco-elastic deviatoric stress (${\mathit{\tau }}_{\mathrm{II}}^{\mathrm{ve}}$) can be obtained by solving the following scalar nonlinear residual equation using the Newton-Raphson method (see Eq. A7):

$$\begin{array}{}\text{(37)}& r^{\mathrm{ve}}={\dot{\mathit{\epsilon }}}_{\mathrm{II}}^{*}-A_{\mathrm{L}}{\mathit{\tau }}_{\mathrm{II}}^{\mathrm{ve}}-A_{\mathrm{N}}{\left({\mathit{\tau }}_{\mathrm{II}}^{\mathrm{ve}}\right)}^n=\mathrm{0},\end{array}$$

where the ${\dot{\mathit{\epsilon }}}_{\mathrm{II}}^{*}$ is the effective deviatoric strain rate that depends on the elastic stresses from the previous time step

$$\begin{array}{}\text{(38)}& {\dot{\mathit{\epsilon }}}_{ij}^{*}={\dot{\mathit{\epsilon }}}_{ij}+{\displaystyle \frac{{\mathit{\tau }}_{ij}^{*}}{\mathrm{2}G\mathrm{\Delta }t}},\end{array}$$

and A_L is the effective linear stress term prefactor, which combines elastic shear and diffusion creep constants:

$$\begin{array}{}\text{(39)}& A_{\mathrm{L}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}G\mathrm{\Delta }t}}+A_{\mathrm{D}}.\end{array}$$

Here ${\mathit{\tau }}_{ij}^{*}$ denotes the rotated deviatoric stress from the previous time step (see the Appendix A for more details). Note that stress rotation terms are usually computed using an alternative incrementally-objective scheme (Thielmann et al., 2015; Gerya, 2019), which produces asymptotically correct results for finite time steps. The 3D generalization of this algorithm was originally provided by Rubinstein and Atluri (1983). The detailed explanation why using the Jaumann objective stress rate does not lead to the spurious stress oscillation for the rocks is provided in Thielmann et al. (2015). The derivative of the visco-elastic scalar residual with respect to the magnitude of the visco-elastic deviatoric stress is readily evaluated as:

$$\begin{array}{}\text{(40)}& {\displaystyle \frac{\partial r^{\mathrm{ve}}}{\partial {\mathit{\tau }}_{\mathrm{II}}^{\mathrm{ve}}}}=-A_{\mathrm{L}}-n{A}_{\mathrm{N}}{\left({\mathit{\tau }}_{\mathrm{II}}^{\mathrm{ve}}\right)}^{n-\mathrm{1}}.\end{array}$$

To ensure robust convergence of the nonlinear iteration we use the following initial guess:

$$\begin{array}{}\text{(41)}& {\mathit{\tau }}_{\mathrm{II}}^{\mathrm{ve}}\approx {\left({\displaystyle \frac{\mathrm{1}}{{\mathit{\tau }}_{\mathrm{L}}}}+{\displaystyle \frac{\mathrm{1}}{{\mathit{\tau }}_{\mathrm{N}}}}\right)}^{-\mathrm{1}},{\mathit{\tau }}_{\mathrm{L}}={\displaystyle \frac{{\dot{\mathit{\epsilon }}}_{\mathrm{II}}^{*}}{A_{\mathrm{L}}}},{\mathit{\tau }}_{\mathrm{N}}={\left({\displaystyle \frac{{\dot{\mathit{\epsilon }}}_{\mathrm{II}}^{*}}{A_{\mathrm{N}}}}\right)}^{\mathrm{1}/n},\end{array}$$

which represents the quasi-harmonic average of two closed-form solutions obtained for each isolated creep mechanism.

During the second stage, the yield function is evaluated using the trial stresses, i.e.: $F=F\left({\mathit{\tau }}_{\mathrm{II}}^{\mathrm{ve}},p^{*}\right)$. If the trial yield function is not violated (F<0), the visco-elastic solution is accepted, i.e.: ${\mathit{\tau }}_{\mathrm{II}}={\mathit{\tau }}_{\mathrm{II}}^{\mathrm{ve}}$, otherwise the viscoplasticity stress constraints are enforced (corrector stage). This step requires solving a coupled system of nonlinear equations with the residual vector (r) incorporating deviatoric (Eq. A7), continuity (Eq. A14) and viscoplastic constraint residuals (Eq. A17). The corresponding solution vector (x) contains the effective deviatoric stress (τ_II), the local pressure (p) and the viscoplastic multiplier ($\dot{\mathit{\lambda }}$), as primary unknowns:

$$\begin{array}{}\text{(42)}& \begin{array}{rl}& \mathit{r}=\left[\begin{array}{rl}& {\dot{\mathit{\epsilon }}}_{\mathrm{II}}^{*}-A_{\mathrm{L}}{\mathit{\tau }}_{\mathrm{II}}-A_{\mathrm{N}}{\left({\mathit{\tau }}_{\mathrm{II}}\right)}^n-\dot{\mathit{\lambda }}{A}_{\mathit{\tau }}\\ & {\displaystyle \frac{p-p^{*}}{K\mathrm{\Delta }t}}-\dot{\mathit{\lambda }}{A}_{\mathrm{p}}\\ & \langle F\rangle -\dot{\mathit{\lambda }}{\mathit{\eta }}^{\mathrm{vp}}\end{array}\right]=\left[\begin{array}{c}\mathrm{0}\\ \mathrm{0}\\ \mathrm{0}\end{array}\right],\\ & \mathit{x}=\left[\begin{array}{c}{\mathit{\tau }}_{\mathrm{II}}\\ p\\ \dot{\mathit{\lambda }}\end{array}\right].\end{array}\end{array}$$

Here the deviatoric and volumetric plastic constants, respectively, are formulated as follows:

$$\begin{array}{}\text{(43)}& \begin{array}{rl}& \left(A_{\mathit{\tau }},A_{\mathrm{p}}\right)=\\ & \left\{\begin{array}{ll}\left(\frac{\mathrm{1}}{\mathrm{2}},k_q\right),& {\mathit{\tau }}_{\mathrm{II}}\left(p_q-p_{\mathrm{d}}\right)\ge {\mathit{\tau }}_{\mathrm{d}}\left(p_q-p\right)\\ \left(\frac{b{\mathit{\tau }}_{\mathrm{II}}}{\mathrm{2}{\widehat{R}}_q},-\frac{b(p-p_q)}{{\widehat{R}}_q}\right),& \mathrm{otherwise}\end{array}\right..\end{array}\end{array}$$

A Newton-Raphson scheme is employed to update the solution until a sufficiently small tolerance is obtained:

$$\begin{array}{}\text{(44)}& {\mathit{x}}_{k+\mathrm{1}}={\mathit{x}}_k-\mathit{\alpha }{\mathbf{J}}_k^{-\mathrm{1}}{\mathit{r}}_k,∥{\mathit{r}}_k∥<\mathrm{tol},\end{array}$$

where k is the iteration index, J the local Jacobian matrix of the local stress update, and ${\mathit{\alpha }}_{\min }<\mathit{\alpha }<\mathrm{1}$ the step length that is optimally selected by a simple back-tracking line-search algorithm guided by the Armijo rule (Armijo, 1966) to ensure that the nonlinear residual is sufficiently reduced between the iterations, i.e:

$$\begin{array}{}\text{(45)}& ∥{\mathit{r}}_{k+\mathrm{1}}∥\le \mathit{\gamma }∥{\mathit{r}}_k∥,\end{array}$$

where $\mathrm{0}<\mathit{\gamma }<\mathrm{1}$ is the residual reduction parameter, which is typically assigned to a rather large value γ=0.9. To ensure progress of the nonlinear iteration, the step length is bounded by a lower threshold α_min. We explicitly emphasize the importance of the line search for the stability of the local iteration algorithms. Despite having a completely smooth yield surface and flow potential, the local iterations can still produce closed loops in the stress space that never lead to convergence when full steps are applied (Fig. 3a). A line search algorithm resolves this issue (Fig. 3c).

Figure 3 Convergence pattern of the nonlinear local iterations for the corrector stage (η^vp=0) when the initial guess is relatively far away from the yield surface. (a) Full Newton steps are applied (α=1). Closed loops and divergence of the algorithm is observed. (b) Damped Newton steps are applied throughout (α=0.1), which results in a smooth pattern, but in a relatively large number of iterations. (c) The step length is controlled by a line search algorithm (γ=0.9, α_min=0.1). Fast convergence is obtained.

To initialize the nonlinear iteration we use the trial visco-elastic stresses and assign the viscoplastic multiplier to zero:

$$\begin{array}{}\text{(46)}& {\mathit{x}}_{\mathrm{0}}=\left[\begin{array}{c}{\mathit{\tau }}_{\mathrm{II}}^{\mathrm{ve}}\\ p^{*}\\ \mathrm{0}\end{array}\right].\end{array}$$

Finally, the Jacobian matrix can be obtained by differentiating the residuals with respect to the unknowns:

$$\begin{array}{}\text{(47)}& \mathbf{J}=\left[\begin{array}{ccc}-A_{\mathrm{L}}-n{A}_{\mathrm{N}}{\left({\mathit{\tau }}_{\mathrm{II}}\right)}^{n-\mathrm{1}}-\dot{\mathit{\lambda }}\frac{\partial A_{\mathit{\tau }}}{\partial {\mathit{\tau }}_{\mathrm{II}}}& -\dot{\mathit{\lambda }}\frac{\partial A_{\mathit{\tau }}}{\partial p}& -A_{\mathit{\tau }}\\ -\dot{\mathit{\lambda }}\frac{\partial A_{\mathrm{p}}}{\partial {\mathit{\tau }}_{\mathrm{II}}}& \frac{\mathrm{1}}{K\mathrm{\Delta }t}-\dot{\mathit{\lambda }}\frac{\partial A_{\mathrm{p}}}{\partial p}& -A_{\mathrm{p}}\\ \frac{\partial F}{\partial {\mathit{\tau }}_{\mathrm{II}}}& \frac{\partial F}{\partial p}& -{\mathit{\eta }}^{\mathrm{vp}}\end{array}\right],\end{array}$$

where the required stress derivatives of the plastic constants (A_τ), (A_p) and the yield function (F) are given by Eqs. (A9), (A16), and (A18), respectively. After convergence of the local iterations, the deviatoric stresses are obtained as follows:

$$\begin{array}{}\text{(48)}& {\mathit{\tau }}_{ij}={\mathit{\tau }}_{\mathrm{II}}{n}_{ij},\end{array}$$

where the normalized deviatoric direction tensor is defined as:

$$\begin{array}{}\text{(49)}& n_{ij}={\displaystyle \frac{{\dot{\mathit{\epsilon }}}_{ij}^{*}}{{\dot{\mathit{\epsilon }}}_{\mathrm{II}}^{*}}}.\end{array}$$

The corresponding update of the accumulated deviatoric viscoplastic strain is given by:

$$\begin{array}{}\text{(50)}& \mathit{\kappa }={\mathit{\kappa }}_n+\mathrm{\Delta }\mathit{\kappa },\mathrm{\Delta }\mathit{\kappa }=\dot{\mathit{\lambda }}{A}_{\mathit{\tau }}\mathrm{\Delta }t.\end{array}$$

To simplify the presentation of the material model in this paper we apply the strain softening explicitly between the time steps. During the local stress update, the yield function is evaluated using the κ values in the beginning of the time step, i.e. F(κ_n). The plastic strength parameters are instantly updated according the Eqs. (50) and (24) once the converged values of Δκ are computed. However, the presented algorithm can be easily modified to include the viscoplastic strain evolution in the local stress update in a coupled manner. This would imply updating Δκ during the iterations and using F(κ) instead. Finally, the accumulated volumetric viscoplastic strain is updated as follows:

$$\begin{array}{}\text{(51)}& \mathit{\chi }={\mathit{\chi }}_n+\mathrm{\Delta }\mathit{\chi },\mathrm{\Delta }\mathit{\chi }=\dot{\mathit{\lambda }}{A}_{\mathrm{p}}\mathrm{\Delta }t.\end{array}$$

3.7 Finite element formulation

We adopt the conforming Crouzeix-Raviart triangular finite element (Crouzeix and Raviart, 1973) to discretize the momentum and continuity residual equations (Eqs. 32 and 34) in the framework of a standard Galerkin procedure. The element is LBB-stable and behaves robustly in practice. Since significant accuracy deterioration was previously reported for the isoparametric finite elements with curvilinear edges (e.g. Lee and Bathe, 1993), we adopt a subparametric formulation and keep the edges of the elements straight. For the accurate evaluation of the element integrals we use an efficient Gaussian quadrature rule specifically designed for triangular elements (Dunavant, 1985). The 2D element has 6 integration points, where all the constitutive equations are evaluated and all the history variables are stored. The details of the finite element discretization are illustrated in Fig. 4.

Figure 4 Conforming triangular Crouzeix-Raviart finite element ($P_{\mathrm{2}}^{+}\times P_{-\mathrm{1}}$). Displacement increments (Δu_e) are interpolated using quadratic shape functions enhanced with a bubble function in the central node. Pressure increments ($\mathbf{\Delta }{\mathit{p}}_e^{*}$) are interpolated using linear shape functions and assumed to be discontinuous between elements. Element coordinates (x_e) are stored only in the element corners, and edges remain straight (subparametric element). History variables and material parameters are stored in the integration points.

In the reminder of this section we summarize the details of the finite element discretization using the standard matrix notation (Zienkiewicz and Taylor, 2000). The primary variables are interpolated in the integration points as follows:

$$\begin{array}{}\text{(52)}& \mathbf{\Delta }\mathit{u}={\mathbf{N}}_u\mathbf{\Delta }{\mathit{u}}_e,\mathrm{\Delta }{p}^{*}={\mathbf{N}}_p\mathbf{\Delta }{\mathit{p}}_e^{*},\end{array}$$

where Δu_e and $\mathbf{\Delta }{\mathit{p}}_e^{*}$ are the displacement and the trial pressure increment vectors of the element, respectively, and the corresponding interpolation matrices are given by:

$$\begin{array}{}\text{(53)}& \begin{array}{rl}& {\mathbf{N}}_u=\left[\begin{array}{ccccc}{N}_{u\mathrm{1}}& \mathrm{0}& \mathrm{\dots }& N_{u\mathrm{7}}& \mathrm{0}\\ \mathrm{0}& N_{u\mathrm{1}}& \mathrm{\dots }& \mathrm{0}& N_{u\mathrm{7}}\end{array}\right],\\ & {\mathbf{N}}_p=\left[\begin{array}{ccc}{N}_{p\mathrm{1}}& \mathrm{\dots }& N_{p\mathrm{3}}\end{array}\right].\end{array}\end{array}$$

Here N_ui and N_pi are respectively the displacement shape function and the pressure shape function of the ith node.

We use Voigt notation to represent the second order symmetric tensors in the vector form. Hence for the total strain increment (Δϵ), the deviatoric strain increment (Δε), and the deviatoric stress (τ) we can write:

$$\begin{array}{}\text{(54)}& \mathbf{\Delta }\mathit{ϵ}=\left[\begin{array}{c}\mathrm{\Delta }{\mathit{ϵ}}_{xx}\\ \mathrm{\Delta }{\mathit{ϵ}}_{yy}\\ \mathrm{\Delta }{\mathit{ϵ}}_{zz}\\ \mathrm{\Delta }{\mathit{\gamma }}_{xy}\end{array}\right],\mathbf{\Delta }\mathit{\epsilon }=\left[\begin{array}{c}\mathrm{\Delta }{\mathit{\epsilon }}_{xx}\\ \mathrm{\Delta }{\mathit{\epsilon }}_{yy}\\ \mathrm{\Delta }{\mathit{\epsilon }}_{zz}\\ \mathrm{\Delta }{\mathit{\epsilon }}_{xy}\end{array}\right],\mathit{\tau }=\left[\begin{array}{c}{\mathit{\tau }}_{xx}\\ {\mathit{\tau }}_{yy}\\ {\mathit{\tau }}_{zz}\\ {\mathit{\tau }}_{xy}\end{array}\right],\end{array}$$

where γ_xy is the engineering shear strain, which is twice larger than the corresponding tensor component, i.e. γ_xy=2 ϵ_xy. The total strain increment vector is computed as:

$$\begin{array}{}\text{(55)}& \mathbf{\Delta }\mathit{ϵ}=\mathbf{B}\mathbf{\Delta }{\mathit{u}}_e,\mathbf{B}=\left[\begin{array}{ccccc}\frac{\partial N_{u\mathrm{1}}}{\partial x}& \mathrm{0}& \mathrm{\dots }& \frac{\partial N_{u\mathrm{7}}}{\partial x}& \mathrm{0}\\ \mathrm{0}& \frac{\partial N_{u\mathrm{1}}}{\partial y}& \mathrm{\dots }& \mathrm{0}& \frac{\partial N_{u\mathrm{7}}}{\partial y}\\ \mathrm{0}& \mathrm{0}& \mathrm{\dots }& \mathrm{0}& \mathrm{0}\\ \frac{\partial N_{u\mathrm{1}}}{\partial y}& \frac{\partial N_{u\mathrm{1}}}{\partial x}& \mathrm{\dots }& \frac{\partial N_{u\mathrm{7}}}{\partial y}& \frac{\partial N_{u\mathrm{7}}}{\partial x}\end{array}\right].\end{array}$$

Here B is the differential operator matrix that contains the derivatives of the displacement shape functions with respect to the global coordinates. Note that the third row explicitly enforces the plane strain kinematical constraint ϵ_zz=0. The deviatoric and volumetric strain increments, respectively, are given by:

$$\begin{array}{}\text{(56)}& \mathbf{\Delta }\mathit{\epsilon }={\mathbf{I}}^D\mathbf{\Delta }\mathit{ϵ},\mathrm{\Delta }\mathit{\theta }={\mathit{m}}^T\mathbf{\Delta }\mathit{ϵ},\end{array}$$

where the corresponding projection matrices are as follows:

$$\begin{array}{}\text{(57)}& {\mathbf{I}}^D=\mathbf{I}-{\displaystyle \frac{\mathrm{1}}{\mathrm{3}}}\mathit{m}{\mathit{m}}^T,\mathbf{I}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left[\begin{array}{cccc}\mathrm{2}& & & \\ & \mathrm{2}& & \\ & & \mathrm{2}& \\ & & & \mathrm{1}\end{array}\right],\mathit{m}=\left[\begin{array}{c}\mathrm{1}\\ \mathrm{1}\\ \mathrm{1}\\ \mathrm{0}\end{array}\right].\end{array}$$

The history variables from the previous time step are assumed to be stored directly in the integration points illustrated as diamonds in Fig. 4. These include the rotated history stress (τ^*), the converged true pressure (p_n), and the accumulated deviatoric (κ_n) and volumetric (χ_n) viscoplastic strains. The following estimate for the effective deviatoric strain rate is readily available:

$$\begin{array}{}\text{(58)}& {\dot{\mathit{\epsilon }}}^{*}={\displaystyle \frac{\mathbf{\Delta }\mathit{\epsilon }}{\mathrm{\Delta }t}}+{\displaystyle \frac{{\mathit{\tau }}^{*}}{\mathrm{2}G\mathrm{\Delta }t}}.\end{array}$$

Combining the Eqs. (36) and (52) we can express the total trial pressure in the integration point as:

$$\begin{array}{}\text{(59)}& p^{*}=p_n+{\mathbf{N}}_p\mathbf{\Delta }{\mathit{p}}_e^{*}.\end{array}$$

Finally, the scalar norm of the effective deviatoric strain rate (${\dot{\mathit{\epsilon }}}_{\mathrm{II}}^{*}$), and the unit deviatoric direction (n) can be evaluated as:

$$\begin{array}{}\text{(60)}& {\dot{\mathit{\epsilon }}}_{\mathrm{II}}^{*}=\sqrt{{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left({\dot{\mathit{\epsilon }}}_{xx}^{*\mathrm{2}}+{\dot{\mathit{\epsilon }}}_{yy}^{*\mathrm{2}}+{\dot{\mathit{\epsilon }}}_{zz}^{*\mathrm{2}}\right)+{\dot{\mathit{\epsilon }}}_{xy}^{*\mathrm{2}}},\mathit{n}={\displaystyle \frac{{\dot{\mathit{\epsilon }}}^{*}}{{\dot{\mathit{\epsilon }}}_{\mathrm{II}}^{*}}}.\end{array}$$

At this stage, we are ready to invoke the local stress update algorithm summarized in Sect. 3.6 to compute the deviatoric stress (τ) and true pressure (p), and hence the total Cauchy stress:

$$\begin{array}{}\text{(61)}& \mathit{\sigma }=\mathit{\tau }-p\mathit{m}.\end{array}$$

The weak forms of the time-discrete momentum and continuity residual equations (Eqs. 32 and 34) are given by:

$$\begin{array}{}\text{(62)}& \begin{array}{rl}{\displaystyle }{\mathit{r}}_u& {\displaystyle }=\underset{V}{\int }{\mathbf{B}}^T\mathit{\sigma }\mathrm{d}V-\underset{V}{\int }{\mathbf{N}}_u^T\mathit{b}\mathrm{d}V-\underset{S}{\int }{\mathbf{N}}_u^T\mathit{t}\mathrm{d}S=\mathrm{0},\\ {\displaystyle }{\mathit{r}}_p& {\displaystyle }=-\underset{V}{\int }{\mathbf{N}}_p^T\mathrm{\Delta }\mathit{\theta }\mathrm{d}V-\underset{V}{\int }{\mathbf{N}}_p^T{\displaystyle \frac{\mathrm{\Delta }{p}^{*}}{K}}\mathrm{d}V=\mathrm{0},\end{array}\end{array}$$

where b=ρg are the gravity forces, and t is the traction vector, which is integrated over the surface. The standard Newton-Raphson iteration is employed to update the primary variables until the global convergence is achieved:

$$\begin{array}{}\text{(63)}& {\left[\begin{array}{c}\mathbf{\Delta }\mathit{u}\\ \mathbf{\Delta }{\mathit{p}}^{*}\end{array}\right]}_{k+\mathrm{1}}={\left[\begin{array}{c}\mathbf{\Delta }\mathit{u}\\ \mathbf{\Delta }{\mathit{p}}^{*}\end{array}\right]}_k-{\left[\begin{array}{ccc}{\mathbf{K}}_{uu}& & {\mathbf{K}}_{up}\\ {\mathbf{K}}_{pu}& & {\mathbf{K}}_{pp}\end{array}\right]}_k^{-\mathrm{1}}{\left[\begin{array}{c}{\mathit{r}}_u\\ {\mathit{r}}_p\end{array}\right]}_k.\end{array}$$

Here k denotes the iteration index, and the blocks of the coupled tangent matrix are evaluated as follows:

$$\begin{array}{}\text{(64)}& \begin{array}{rlrl}{\displaystyle }{\mathbf{K}}_{uu}& {\displaystyle }=\underset{V}{\int }{\mathbf{B}}^T\mathbf{D}\mathbf{B}\mathrm{d}V,& {\displaystyle }{\mathbf{K}}_{up}& {\displaystyle }=-\underset{V}{\int }{\mathbf{B}}^T\mathit{q}{\mathbf{N}}_p\mathrm{d}V,\\ {\displaystyle }{\mathbf{K}}_{pu}& {\displaystyle }=-\underset{V}{\int }{\mathbf{N}}_p^T{\mathit{m}}^T\mathbf{B}\mathrm{d}V,& {\displaystyle }{\mathbf{K}}_{pp}& {\displaystyle }=-\underset{V}{\int }{\mathbf{N}}_p^T{\displaystyle \frac{\mathrm{1}}{K}}{\mathbf{N}}_p\mathrm{d}V,\end{array}\end{array}$$

where the tangent operators are given by:

$$\begin{array}{}\text{(65)}& \begin{array}{rl}& \mathbf{D}={\displaystyle \frac{\partial \mathit{\sigma }}{\partial \mathbf{\Delta }\mathit{ϵ}}}=\left(\mathrm{2}{\mathit{\eta }}_{\mathrm{eff}}{\mathbf{I}}^D+{\mathit{\beta }}_{\mathrm{1}}\mathit{n}{\mathit{n}}^T+{\mathit{\beta }}_{\mathrm{2}}\mathit{m}{\mathit{n}}^T\right)/\mathrm{\Delta }t,\\ & \mathit{q}={\displaystyle \frac{\partial \mathit{\sigma }}{\partial p^{*}}}={\mathit{\beta }}_{\mathrm{3}}\mathit{n}+{\mathit{\beta }}_{\mathrm{4}}\mathit{m}.\end{array}\end{array}$$

Details of the derivations along with expressions for the coefficients η_eff, and β₁−β₄ are provided in Appendix B. It should be noted that only completely pressure-independent constitutive models (such as von Mises plasticity or nonlinear visco-elasticity) produce a symmetric tangent matrix. In these cases, the stiffness coefficients assume the following values: ${\mathit{\beta }}_{\mathrm{2}}={\mathit{\beta }}_{\mathrm{3}}=\mathrm{0},{\mathit{\beta }}_{\mathrm{4}}=\mathrm{1}$. Any form of pressure dependence in the constitutive model, including the fully-associated plasticity case, introduces non-symmetric terms in the tangent matrix. This property of the mixed two-field formulation contrasts with the classical strain-driven approach which preserves symmetry for the associated pressure-dependent plasticity models.

The global Newton iteration is terminated as soon as the following stopping criterion is satisfied:

$$\begin{array}{}\text{(66)}& {\displaystyle \frac{\left|\left|{\mathit{f}}^{\mathrm{int}}-{\mathit{f}}^{\mathrm{ext}}\right|\right|}{\left|\left|{\mathit{f}}^{\mathrm{ext}}\right|\right|}}<\mathrm{tol},\end{array}$$

where the partition of the momentum residual vector into the external and internal forces is defined as:

$$\begin{array}{}\text{(67)}& {\mathit{f}}^{\mathrm{int}}=\underset{V}{\int }{\mathbf{B}}^T\mathit{\sigma }\mathrm{d}V,{\mathit{f}}^{\mathrm{ext}}=\underset{V}{\int }{\mathbf{N}}_u^T\mathit{b}\mathrm{d}V+\underset{S}{\int }{\mathbf{N}}_u^T\mathit{t}\mathrm{d}S=\mathrm{0}\end{array}$$

There is no need to explicitly check the residual of the continuity equation, since it is typically satisfied to a very low tolerance at every global iteration.

After achieving global convergence, velocities, deviatoric strain rates, and volumetric strain rates, respectively, can be computed for postprocessing purposes using:

$$\begin{array}{}\text{(68)}& \mathit{v}={\displaystyle \frac{\mathbf{\Delta }\mathit{u}}{\mathrm{\Delta }t}},\dot{\mathit{\epsilon }}={\displaystyle \frac{\mathbf{\Delta }\mathit{\epsilon }}{\mathrm{\Delta }t}},\dot{\mathit{\theta }}={\displaystyle \frac{\mathrm{\Delta }\mathit{\theta }}{\mathrm{\Delta }t}}.\end{array}$$

3.8 Software implementation

To test the applicability and robustness of the yield surface developed in this work, we implemented the algorithms presented in Sect. 3.6 and 3.7 in a compact Python finite element code (GeoTech2D) that employs LBB-stable triangular Crouziex-Raviart elements. We use an efficient array implementation from the NumPy package (Harris et al., 2020). The finite element matrices are initially stored in coordinate format and subsequently assembled in a compressed sparse column format using the SciPy package (Virtanen et al., 2020). Dabrowski et al. (2008) implemented a similar approach and demonstrated that it has superior efficiency compared to a direct assembly. The default sparse direct solver provided by the SciPy package is used to solve the linearized system of equations. We do not utilize the block structure of the tangent matrix by adopting efficient approaches such as Powell-Hestenes iterations (e.g. Dabrowski et al., 2008). Since we only consider elastically compressible cases the robustness of the direct solver is sufficient to perform a coupled factorization. To improve the robustness and speed of the code, iterative solvers could be added at later stage. Refer to the code availability section to obtain the instructions on how to access and run the code.

To generate the unstructured triangular grids we use the Triangle mesh generator (Shewchuk, 1996), which we invoke via the Python interface implemented in the MeshPy package (Kloeckner et al., 2022). In the code availability section we provide the links to access all the scripts that we used to generate the grids in this paper along with an example script that demonstrates the basic features of the Triangle interface in MeshPy, such as material regions, holes, and boundary markers. We use the pyEVTK package (Herrera, 2021) to store the simulation results in the Visualization Toolkit (VTK) format (Schroeder et al., 2006), which can be visualized with the ParaView package (Ahrens et al., 2005). Whenever necessary we use the scientific color maps (Crameri et al., 2020), which prevent visual distortion of the data.

4 Benchmarks and applications

In this section we test the application of the developed plasticity model to a number of relevant problems that involve both mode-I and mode-II plastic failure. The material parameters, temporal and spatial discretization details, domain geometry, and boundary conditions are described in Table 1. Every problem we consider here is summarized in a single column of this table, which is labeled with a keyword concisely describing the essence of the setup (e.g. Crust or Tensile, etc). References to the figures that demonstrate basic results of a corresponding numerical simulation are also provided. Parameters related to different scenarios of the considered problems (e.g. extension vs. compression) are separated by a slash.

Table 1 Material, discretization, and boundary condition parameters employed in presented models.

Note: L – domain length, H – domain height, A_△ – target triangular element area, ${\dot{\mathit{ϵ}}}_{ij}^{\mathrm{bg}}$ – background strain rates.

4.1 0D stress integration test

The correctness of the stress integration algorithm implementation can be demonstrated by performing the so-called 0D deformation experiment (integration point test). The major purpose of this test is to show that local iterations adequately calculate stresses during the switch from the visco-elastic deformation into mode-I failure, as well as during the transition between mode-I and mode-II. Material properties are assumed to be homogeneous and magnitudes of the background strain rates components are selected such that three scenarios are generated, namely: (i) constant volumetric strain rate is applied $\dot{\mathit{\theta }}\ne \mathrm{0}$ (volumetric extension test) (Fig. 5a), (ii) constant deviatoric strain rate is applied ${\dot{\mathit{\epsilon }}}_{\text{II}}\ne \mathrm{0}$ (deviatoric shear test) (Fig. 5b), and (iii) combination of both strain rates is applied (mixed strain test) (Fig. 5c). For simplicity we consider perfect plasticity case, i.e. we ignore softening and deactivate viscoplastic regularization. All three scenarios are initialized with zero stresses and subsequently, either pressure or deviatoric stress, or both, evolve until mode-I failure is activated. At later stages the stress evolution only occurs along the yield surface. For the constant volumetric strain rate case, the pressure simply stops changing when it reaches the tensile strength of the material (Fig. 5a, d). In the deviatoric and mixed cases the mode-I failure eventually switches to the mode-II when stress evolution path passes through the delimiter point between the yield surface segments (Fig. 5b, c, d). Due to the nonzero dilatation angle, the pressure continues to grow even after switching to the mode-II regime. This behavior is caused by a combination of the nonzero volumetric plastic strain and overall incompressibility constraint ($\dot{\mathit{\theta }}=\mathrm{0}$), which results in elastic compression and associated pressure increase. In general the stress integration tests deliver the expected results for the considered loading scenarios.

Figure 5 Temporal evolution of the effective deviatoric stress and pressure for the 0D stress integration tests. (a) Volumetric extension test $\dot{\mathit{\theta }}=\mathrm{7}\times {\mathrm{10}}^{-\mathrm{15}}{\mathrm{s}}^{-\mathrm{1}}$ (b) Deviatoric shear test ${\dot{\mathit{\epsilon }}}_{\text{II}}=\mathrm{7}\times {\mathrm{10}}^{-\mathrm{14}}{\mathrm{s}}^{-\mathrm{1}}$. (c) Mixed strain test (a combination of the volumetric extension and deviatoric shear used in a and b). (d) Stress evolution patterns in the meridional profile for all three tests. Both mode-I and transition between mode-I and mode-II are tested.

4.2 Viscoplastic regularization test

In the next example, we test the effect of viscoplastic regularization in a 2D setup that develops localized zones by either mode-I or mode-II plasticity. Two combinations of background strain rate loading is considered that correspond either to the uniaxial restrained extension (Fig. 6a, b), or to the pure shear (Fig. 6c, d). For the extension case we plot the accumulated volumetric viscoplastic strain (χ), which indicate loci of the tensile failure zones, whereas the deviatoric counterpart (κ) is plotted for the pure share case, which corresponds to the shear bands. The test cases are performed for different mesh resolutions. As expected, the non-regularized plasticity models demonstrates severe mesh dependency of both mode-I and mode-II plasticity cases, which results in localization zones that are close to one element in width and that have a maximum strain that depends on the mesh resolution (Fig. 6a, c). With a suitably tuned viscoplastic regularization viscosity, on the other hand, mesh independent results are obtained (Fig. 6b, d), similar to earlier results for mode-II plasticity (Duretz et al., 2019, 2023). It should be noted that the Perzyna-type regularization we employ affects both deviatoric and volumetric components of viscoplastic strain rate via the viscoplastic multiplier ($\dot{\mathit{\lambda }}$) (Eqs. 5, 9, and 11), and hence acts on mode-I and mode-II plasticity simultaneously. As previously reported by Duretz et al. (2019), one of the apparent advantages of the regularized models is its greatly improved convergence rate compared to the non-regularized models. In particular a few high resolution non-regularized models we tested here simply failed to converge due to the residual stagnation (Fig. 6a, c). In contrast to that, all the regularized models converged successfully within a prescribed number of iterations for both mode-I and mode-II cases (Fig. 6b, d). The typical convergence profiles of the nonlinear iterations are discussed in more detail in Sect. 4.6.

Figure 6 Results of viscoplastic regularization tests. (a) Tensile localization test without regularization after 10 kyr. (b) Tensile localization test with regularization (η^vp=10¹⁹) after 100 kyr. (c) Shear localization test without regularization after 1 kyr. (d) Shear localization test with regularization (${\mathit{\eta }}^{\mathrm{vp}}=\mathrm{5}\times {\mathrm{10}}^{\mathrm{18}}$) after 3 kyr. Curves on the cross-sections are labeled with target triangular element areas (mesh resolutions). Spatial distribution plots of χ and κ correspond to the highest resolution tested for each case.

4.3 Strain localization in the brittle crust

Next, we consider the deformation of the upper crust in either extensional or compressional setting with an elasto-viscoplastic rheology. This problem represents a standard benchmark used in the geodynamics community to test various plasticity model implementations (see e.g. Kaus, 2010, and references therein). To facilitate the localization we use strain softening for cohesion parameter which is initiated by random accumulated deviatoric viscoplastic strain perturbations that are clustered around the central upper part of the domain. Subsequent deformation manifests itself in formation of multiple incipient localization zones (Fig. 7a), among which only few develop into fault zones at later stages (Fig. 7b). Under extension, mode-I plastic failure occurs close to the free surface, which is visible as vertically oriented high strain rate zones in Fig. 7a. With increasing depth and confining stress the localization switches again to mode-II plasticity. Under compression, the mode-I failure does not occur as mean stresses never become extensional. In general, it takes more time and horizontal strain to develop the localization zones in the compressional setting. In the considered setup the formation of the faults finishes at about 2.4 kyr under extension (Fig. 7b), whereas it still continues after about 7 kyr under compression (Fig. 7c). This can be explained by the increased compressive strength due to the increased dynamic pressure caused by the tectonic shortening. The orientation angles of the shear bands produced by the numerical models is the major observation parameter that is compared against the theoretical estimates for this benchmark problem. We refer to Kaus (2010) for the overview of the existing estimates as well as for the extensive discussion regarding the parameters that influence these orientations. Typically the result produced by the numerical models follow the theoretical estimates of Arthur et al. (1978), given by:

$$\begin{array}{}\text{(69)}& {\mathit{\alpha }}_{\mathrm{S}}={\displaystyle \frac{\mathit{\pi }}{\mathrm{4}}}\pm {\displaystyle \frac{\mathit{\phi }+\mathit{\psi }}{\mathrm{4}}}.\end{array}$$

Here α_S stands for the dip angles of the shear bands. The positive sign in this expression corresponds to the extensional setup, whereas the negative sign is attributed to compression. In general, normal faults are predicted to localize much steeper compared to the thrust faults. Detailed inspection of Fig. 7 reveals that Arthur's estimates are reproduced quite accurately in the presented simulations.

Figure 7 Spatial distribution of the deviatoric strain rate for the crustal strain localization problem. (a) Extension case after 1.8 kyr (incipient localization) and 2.4 kyr (localized faults). (b) Compression case after 7 kyr. Approximate shear band orientation angles are explicitly indicated. Note that the top of the domain is stress-free (representing the Earth's surface).

4.4 Tensile failure zone propagation

In this experiment, we consider propagation of a localized tensile failure zone induced by elevated fluid pressure. According to the theory of poroelasticity (Biot, 1941) the pressure variable in the momentum equation (Eq. 10) must be replaced with the total pressure given by:

$$\begin{array}{}\text{(70)}& p_{\mathrm{total}}=p+{\mathit{\alpha }}_{\mathrm{B}}{p}_{\mathrm{f}},\end{array}$$

where p_f denotes the fluid pressure in the rock pores, α_B is the Biot-Willis coefficient, which depends on the bulk modulii of the rock matrix and solid grains (Biot and Willis, 1957), and p is interpreted as the effective pressure that controls the deformation and failure of the porous rock. Apart from this single modification, the rest of the theory and notation in this paper remains unaltered. For the current test we set the Biot-Willis coefficient to one, which simply corresponds to Terzaghi's effective stress principle (Terzaghi, 1925). For all other tests in this paper, we set the Biot-Willis coefficient to zero and thus ignore the effect of fluid pressure.

The initial fluid pressure distribution is assumed to be homogeneous. We override it with a local perturbation of 180 MPa at the center point of the bottom boundary. The perturbation is assumed to decay around the center point following a Gaussian distribution with standard deviation of 200 m in the horizontal direction and 500 m in the vertical direction. The initial fluid pressure perturbation is sufficiently large to initiate tensile failure. Subsequently the excessive fluid pressure is propagated to the regions where the accumulated volumetric viscoplastic strain exceeds the prescribed limit, i.e.: χ>χ_lim, which is assigned to a relatively low value e.g. 10⁻⁴. This explicit procedure effectively mimics Darcy flow of the fluid into the enhanced permeability zone induced by hydraulic fracturing, even when we do not model the actual underlying physical processes related to fluid flow, in order to focus on the mechanical aspects of tensile failure.

The positive feedback between the plastic strain and the fluid pressure leads to a gradual propagation of the failure zone in the vertical direction (Fig. 8a). The propagation velocity of the zone is increasing as the localization progresses towards the surface. The fluid pressure in the zone remains almost the same as in the original perturbation due to the lack of viscous dissipation (Fig. 8b), while the total pressure rapidly decreases due to decreasing overburden weight at shallower depths. Hence the effective pressure becomes tensile and grows in magnitude, causing faster rock failure at the tip of the zone and therefore faster propagation rates. The velocities at the boundaries of the domain are essentially controlled by the imposed background strain rates. However, near the tip of the zone they reach values that are orders of magnitude larger than at the boundaries (Fig. 8c), suggesting that the dynamics of the system is mainly driven by the brittle failure and associated elastic unloading. The animation of the time evolution of the horizontal velocity field is provided in the video supplement.

Figure 8 Results of tensile failure zone propagation problem. (a) Accumulated volumetric viscoplastic strain distribution after 12 years (onset) and 15 years (breakthrough). (b) Fluid pressure distribution after 12 and 15 years. (c) Horizontal velocity distribution after 12 years.

It should be noted that this setup is not intended to be interpreted as dyke propagation through the brittle crust even though it resembles it rather closely. The reasons for this are numerous. We ignore thermal effects, do not represent realistic magma viscosities and densities in the model, describe fluid flow by a simplified parametrization, etc. In reality, dikes and hydraulic fractures form on much smaller length- and time-scales, compared to what is represented here. More realistic dyke propagation models would require much finer numerical resolution.

4.5 Brittle-ductile transition

The tests so far focused on elasto-viscoplastic rheologies. Yet, our numerical formulation can also deal with linear and nonlinear creep laws combined with elasto-viscoplastic failure models. To demonstrate that, our final test considers a similar crustal extension setup we already introduced in Sect. 4.3. However this time we significantly increase the horizontal and vertical spans of the domain to be 100 × 25 km. The domain is discretized with a variable grid resolution. The central upper part of the domain with a size of 40 × 7 km has a target element area of 3600 m² (60 m average mesh size), which is just slightly greater than in the crustal extension setup. The background element area is assigned to 25 000 m² (500 m average mesh size). Within a 20 km transition zone the element size is assumed to linearly grow from the refined value to a background value to prevent sharp contrast in grid resolution. We assigned a background geothermal gradient to 20 °C km⁻¹ and extend the crust in the horizontal direction with a constant background strain rate. A combination of the selected background geothermal gradient, a constant extension rate and a temperature- and stress-dependent dry upper crust quartzite rheology (Schmalholz et al., 2009) ensures that the brittle-ductile transition occurs within the model domain, which is confirmed by a typical strength envelope of the crust (Fig. 9c). In the brittle part of the crust we again obtain a combination of mode-I and mode-II plastic failure zones, which develop in the central part of the domain (Fig. 9a, b). The localization is triggered by similar random perturbations of the accumulated deviatoric viscoplastic strain as we used in the crustal setup (see Sect. 4.3 for the details). Peculiar features that we obtain in our models are the vertical mode-I failure zones that occur close to the free surface where compressive strength is the smallest due to limited overburden stress (Fig. 9b). It should be noted that these features are not reproducible by the models that do not include a tensile yield surface (Popov and Sobolev, 2008; Kaus, 2010; Duretz et al., 2021; Jacquey and Cacace, 2020a, b). The animation of the time evolution of the accumulated deviatoric viscoplastic strain field is provided in the video supplement.

Figure 9 Results of brittle-ductile transition problem. (a) Accumulated deviatoric viscoplastic strain distribution after 33.5 kyr. (b) Accumulated volumetric viscoplastic strain distribution in a near-surface region. (c) Depth distribution of the effective deviatoric stress magnitude in the far field. The top of the domain is stress-free.

4.6 Convergence of the global iterations

Global Newton-Raphson iteration, which we use in this paper, is a powerful technique to solve nonlinear systems of equations. Nevertheless, it is well known that non-regularized plasticity models usually demonstrate very poor convergence rates, residual stagnation, or even divergence caused by the unlimited residual growth. As previously reported by e.g. Duretz et al. (2019), the positive side effect of the viscoplastic regularization is its significantly more robust convergence properties compared to the non-regularized cases. Here, we generally confirm this observation. For a typical crustal extension problem considered in Sect. 4.3 only few (generally less than 10, and normally about 4–5) nonlinear iterations are necessary to achieve a relative tolerance of 10⁻⁵ (Fig. 10). In cases where an overwhelming amount of loading is applied during a single time step, which typically happens when a lot of fault zones are active simultaneously, the residual starts to stagnate and the nonlinear solve fails to converge within a prescribed number of iterations. In this situation, we use an adaptive time stepping, which simply halves the timestep and restarts the nonlinear solution procedure. After a prescribed number of successfully converged reduced times steps, the algorithm attempts to double the time step magnitude again to avoid excessively refined time discretization. It should finally be noted that the robustness of the nonlinear solution can be potentially further improved by a line search technique, similar to what we apply for the local stress update. In this paper, however, we did not explore this possibility for the global equilibrium iterations.

Figure 10 Snapshot of the typical nonlinear global convergence pattern for the crustal extension model (Sect. 4.3). Shown are the logarithm of the relative momentum residual norm vs. the cumulative iteration number. The partition of the momentum residual into internal and external forces is given by Eq. (67). Green dashed line represents the relative tolerance value used to terminate the iterations. Stars indicate converged residual norms of the time steps.

5 Conclusions

We proposed a relatively simple plasticity model that combines a linearized Drucker-Prager shear failure envelope with a circular tensile cap function in a way that ensures a globally smooth transition between these components. Similar models were elaborated before, however, here we additionally eliminate all singularities and discontinuities in the entire stress space for both the composite yield surface, and its corresponding flow potential. This important addition makes our plasticity model a suitable candidate for implementation in the context of Perzyna-type viscoplasticity, which is a key component to obtain regularized solutions that are independent of the mesh resolution. The proposed plasticity model is characterized by a unique orientation of the flow potential gradient at any stress point, which in turn allows to formulate a relatively simple and robust local stress integration algorithm. Nevertheless, we found it necessary to use a line search algorithm to increase the stability of the local iterations and to achieve a practically acceptable level of robustness. Our rheological model also supports an arbitrary combination of viscoplasticity with linear and nonlinear viscous creep mechanisms, elasticity and strain softening. The amount of dilatation in the shear failure envelope does not affect the formulation of the tensile cap surface, which implies that zero dilatation angle cases are explicitly supported.

We implemented the proposed plasticity model in a new unstructured finite element code that uses LBB-stable conforming triangular Crouzier-Raviart elements to treat the near incompressibility constraint. Our mixed two-field formulation employs increments of displacements and pressure as primary variables. We found that treating the globally discretized pressure variable as a trial visco-elastic pressure (trial pressure scheme) performs superior compared to treating it as a true spherical part of the Cauchy stress tensor (true pressure scheme). We also found that displacement increments are advantageous compared to velocities as they help to control kinematical loading by changing the time step magnitude, while posing no restrictions on the type of rheology employed.

A number of 0D and 2D examples are provided that demonstrate the applicability of the proposed plasticity model to a set of practically relevant cases that involve both mode-I and mode-II plasticity. These range from fluid-induced tensile failure zone propagation to strain localization in the visco-elasto-plastic crust. In general, it can be concluded that the algorithms and models presented in this work are quite robust and can be easily implemented in the framework of mixed finite element or staggered-grid finite difference formulations to simulate a wide range of geomechanical applications.

Appendix A: Stress update

Applying the time discretization to the Jaumann stress rate (Eq. 6) we obtain:

$$\begin{array}{}\text{(A1)}& {\stackrel{\diamond }{\mathit{\tau }}}_{ij}={\displaystyle \frac{{\mathit{\tau }}_{ij}-{\left({\mathit{\tau }}_{ij}\right)}_n}{\mathrm{\Delta }t}}+{\left({\mathit{\tau }}_{ik}\right)}_n{\mathit{\omega }}_{kj}-{\mathit{\omega }}_{ik}{\left({\mathit{\tau }}_{kj}\right)}_n.\end{array}$$

Taking the deviatoric projection of plastic flow potential gradient (Eq. 21) gives:

$$\begin{array}{}\text{(A2)}& \mathrm{dev}\left({\displaystyle \frac{\partial Q}{\partial {\mathit{\sigma }}_{ij}}}\right)=B_{\mathit{\tau }}{\mathit{\tau }}_{ij}.\end{array}$$

Upon substitution of both results in the deviatoric constitutive equation (Eq. 5) we finally get:

$$\begin{array}{}\text{(A3)}& \begin{array}{rl}{\dot{\mathit{\epsilon }}}_{ij}-& {\displaystyle \frac{\mathrm{1}}{\mathrm{2}G}}\left[{\displaystyle \frac{{\mathit{\tau }}_{ij}-{\left({\mathit{\tau }}_{ij}\right)}_n}{\mathrm{\Delta }t}}+{\left({\mathit{\tau }}_{ik}\right)}_n{\mathit{\omega }}_{kj}-{\mathit{\omega }}_{ik}{\left({\mathit{\tau }}_{kj}\right)}_n\right]\\ & -A_{\mathrm{D}}{\mathit{\tau }}_{ij}-A_{\mathrm{N}}{\left({\mathit{\tau }}_{\mathrm{II}}\right)}^{n-\mathrm{1}}{\mathit{\tau }}_{ij}-\dot{\mathit{\lambda }}{B}_{\mathit{\tau }}{\mathit{\tau }}_{ij}=\mathrm{0}.\end{array}\end{array}$$

Rearranging the results yields:

$$\begin{array}{}\text{(A4)}& {\dot{\mathit{\epsilon }}}_{ij}^{*}-A_{\mathrm{L}}{\mathit{\tau }}_{ij}-A_{\mathrm{N}}{\left({\mathit{\tau }}_{\mathrm{II}}\right)}^{n-\mathrm{1}}{\mathit{\tau }}_{ij}-\dot{\mathit{\lambda }}{B}_{\mathit{\tau }}{\mathit{\tau }}_{ij}=\mathrm{0},\end{array}$$

where the effective linear visco-elastic creep prefactor is defined as:

$$\begin{array}{}\text{(A5)}& A_{\mathrm{L}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}G\mathrm{\Delta }t}}+A_{\mathrm{D}},\end{array}$$

and the effective deviatoric strain rate and rotated deviatoric stresses from the previous time step, respectively, are given by:

$$\begin{array}{}\text{(A6)}& \begin{array}{rl}& {\dot{\mathit{\epsilon }}}_{ij}^{*}={\dot{\mathit{\epsilon }}}_{ij}+{\displaystyle \frac{{\mathit{\tau }}_{ij}^{*}}{\mathrm{2}G\mathrm{\Delta }t}},\\ & {\mathit{\tau }}_{ij}^{*}={\left({\mathit{\tau }}_{ij}\right)}_n+\mathrm{\Delta }t\left[{\mathit{\omega }}_{ik}{\left({\mathit{\tau }}_{kj}\right)}_n-{\left({\mathit{\tau }}_{ik}\right)}_n{\mathit{\omega }}_{kj}\right].\end{array}\end{array}$$

From Eq. (A4) it becomes obvious that the updated deviatoric stress and the effective strain rate are proportional to each other, which implies that all terms can be replaced with their corresponding scalar norms:

$$\begin{array}{}\text{(A7)}& {\dot{\mathit{\epsilon }}}_{\mathrm{II}}^{*}-A_{\mathrm{L}}{\mathit{\tau }}_{\mathrm{II}}-A_{\mathrm{N}}{\left({\mathit{\tau }}_{\mathrm{II}}\right)}^n-\dot{\mathit{\lambda }}{A}_{\mathit{\tau }}=\mathrm{0}\end{array}$$

where the deviatoric plastic constant is computed as:

$$\begin{array}{}\text{(A8)}& A_{\mathit{\tau }}=B_{\mathit{\tau }}{\mathit{\tau }}_{\mathrm{II}}=\left\{\begin{array}{ll}\frac{\mathrm{1}}{\mathrm{2}},& {\mathit{\tau }}_{\mathrm{II}}\left(p_q-p_{\mathrm{d}}\right)\ge {\mathit{\tau }}_{\mathrm{d}}\left(p_q-p\right)\\ \frac{b{\mathit{\tau }}_{\mathrm{II}}}{\mathrm{2}{\widehat{R}}_q},& \mathrm{otherwise}\end{array}\right..\end{array}$$

The stress derivatives of the deviatoric plastic constant are given by:

$$\begin{array}{}\text{(A9)}& \begin{array}{rl}& \left({\displaystyle \frac{\partial A_{\mathit{\tau }}}{\partial {\mathit{\tau }}_{\mathrm{II}}}},{\displaystyle \frac{\partial A_{\mathit{\tau }}}{\partial p}}\right)=\\ & \left\{\begin{array}{ll}\left(\mathrm{0},\mathrm{0}\right),& {\mathit{\tau }}_{\mathrm{II}}\left(p_q-p_{\mathrm{d}}\right)\ge {\mathit{\tau }}_{\mathrm{d}}\left(p_q-p\right)\\ \left(\frac{b\left({\widehat{R}}_q^{\mathrm{2}}-{\mathit{\tau }}_{\mathrm{II}}^{\mathrm{2}}\right)}{\mathrm{2}{\widehat{R}}_q^{\mathrm{3}}},-\frac{b{\mathit{\tau }}_{\mathrm{II}}(p-p_q)}{\mathrm{2}{\widehat{R}}_q^{\mathrm{3}}}\right),& \mathrm{otherwise}\end{array}\right..\end{array}\end{array}$$

For the continuity equation, we follow essentially the same steps as for the deviatoric equation. By discretizing the pressure time derivative we get:

$$\begin{array}{}\text{(A10)}& {\displaystyle \frac{\mathrm{D}p}{\mathrm{D}t}}{\displaystyle }={\displaystyle \frac{p-p_n}{\mathrm{\Delta }t}}.\end{array}$$

The spherical part of the plastic flow potential gradient (Eq. 21) reads:

$$\begin{array}{}\text{(A11)}& {\displaystyle }\mathrm{tr}\left({\displaystyle \frac{\partial Q}{\partial {\mathit{\sigma }}_{ij}}}\right){\displaystyle }=\mathrm{3}{B}_{\mathrm{p}}.\end{array}$$

Substituting both terms in the continuity equation (Eq. 9) gives:

$$\begin{array}{}\text{(A12)}& \dot{\mathit{\theta }}=-{\displaystyle \frac{p-p_n}{K\mathrm{\Delta }t}}+\mathrm{3}\dot{\mathit{\lambda }}{B}_{\mathrm{p}}.\end{array}$$

Assuming zero viscoplastic multiplier, the volumetric strain rate can be expressed in terms of the trial pressure as:

$$\begin{array}{}\text{(A13)}& \dot{\mathit{\theta }}=-{\displaystyle \frac{p^{*}-p_n}{K\mathrm{\Delta }t}}.\end{array}$$

Combining the right hand sides of both equations and rearranging the result, we finally obtain the local continuity residual:

$$\begin{array}{}\text{(A14)}& {\displaystyle \frac{p-p^{*}}{K\mathrm{\Delta }t}}-\dot{\mathit{\lambda }}{A}_{\mathrm{p}}=\mathrm{0},\end{array}$$

where the volumetric plastic constant is computed as:

$$\begin{array}{}\text{(A15)}& A_{\mathrm{p}}=\mathrm{3}{B}_{\mathrm{p}}=\left\{\begin{array}{ll}{k}_q,& {\mathit{\tau }}_{\mathrm{II}}\left(p_q-p_{\mathrm{d}}\right)\ge {\mathit{\tau }}_{\mathrm{d}}\left(p_q-p\right)\\ -\frac{b(p-p_q)}{{\widehat{R}}_q},& \mathrm{otherwise}\end{array}\right..\end{array}$$

The stress derivatives of the volumetric plastic constant are given by:

$$\begin{array}{}\text{(A16)}& \begin{array}{rl}& \left({\displaystyle \frac{\partial A_{\mathrm{p}}}{\partial {\mathit{\tau }}_{\mathrm{II}}}},{\displaystyle \frac{\partial A_{\mathrm{p}}}{\partial p}}\right)=\\ & \left\{\begin{array}{ll}\left(\mathrm{0},\mathrm{0}\right),& {\mathit{\tau }}_{\mathrm{II}}\left(p_q-p_{\mathrm{d}}\right)\ge {\mathit{\tau }}_{\mathrm{d}}\left(p_q-p\right)\\ \left(\frac{b{\mathit{\tau }}_{\mathrm{II}}(p-p_q)}{{\widehat{R}}_q^{\mathrm{3}}},-\frac{b\left({\widehat{R}}_q^{\mathrm{2}}-(p-p_q{)}^{\mathrm{2}}\right)}{{\widehat{R}}_q^{\mathrm{3}}}\right),& \mathrm{otherwise}\end{array}\right..\end{array}\end{array}$$

Finally, the Perzyna viscoplastic constitutive equation (Eq. 11) can be recast into a residual form:

$$\begin{array}{}\text{(A17)}& \dot{\mathit{\lambda }}={\displaystyle \frac{\langle F\rangle }{{\mathit{\eta }}^{\mathrm{vp}}}}\to \langle F\rangle -\dot{\mathit{\lambda }}{\mathit{\eta }}^{\mathrm{vp}}=\mathrm{0}.\end{array}$$

The stress derivatives of the composite yield surface are given by:

$$\begin{array}{}\text{(A18)}& \begin{array}{rl}& \left({\displaystyle \frac{\partial F}{\partial {\mathit{\tau }}_{\mathrm{II}}}},{\displaystyle \frac{\partial F}{\partial p}}\right)=\\ & \left\{\begin{array}{ll}\left(\mathrm{1},-k\right),& {\mathit{\tau }}_{\mathrm{II}}\left(p_{\mathrm{y}}-p_{\mathrm{d}}\right)\ge {\mathit{\tau }}_{\mathrm{d}}\left(p_{\mathrm{y}}-p\right)\\ \left(\frac{a{\mathit{\tau }}_{\mathrm{II}}}{{\widehat{R}}_y},\frac{a(p-p_{\mathrm{y}})}{{\widehat{R}}_y}\right),& \mathrm{otherwise}\end{array}\right..\end{array}\end{array}$$

Appendix B: Tangent operator

To simplify the notation, we omit the star superscript indicating the effective deviatoric strain rates $\left({\dot{\mathit{\epsilon }}}_{ij}^{*}\right)$ in the following derivations. The linearization of the stress rotation terms, in general, depends on the employed advection algorithm. It can be relatively simple in case that rotation is performed in the integration points or control volumes. However, linearization of the rotation terms performed on the material particles during the advection can be notoriously difficult. Since no general algorithm can be provided we assume that stress rotation terms do not contribute to the linearization, which is a fair assumption for geoscientific applications as rocks reach their ultimate yield stress and break when shear stresses do not exceed about 10 % of the elastic shear modulus. At this stress level, stress rotation terms do not play a crucial role yet (see e.g., Thielmann et al., 2015). We therefore limit the derivations to the linearization of the material nonlinearity.

Directly differentiating the Newton update Eq. (44) with respect to the effective deviatoric strain rate (${\dot{\mathit{\epsilon }}}_{\mathrm{II}}$), and the trial pressure (p^*), assuming a unit step length (α=1), yields the following matrix:

$$\begin{array}{}\text{(B1)}& \mathbf{A}=\left[\begin{array}{cc}\frac{\partial {\mathit{\tau }}_{\mathrm{II}}}{\partial {\dot{\mathit{\epsilon }}}_{\mathrm{II}}}& \frac{\partial {\mathit{\tau }}_{\mathrm{II}}}{\partial p^{*}}\\ \frac{\partial p}{\partial {\dot{\mathit{\epsilon }}}_{\mathrm{II}}}& \frac{\partial p}{\partial p^{*}}\\ \frac{\partial \dot{\mathit{\lambda }}}{\partial {\dot{\mathit{\epsilon }}}_{\mathrm{II}}}& \frac{\partial \dot{\mathit{\lambda }}}{\partial p^{*}}\end{array}\right]=-{\mathbf{J}}^{-\mathrm{1}}\left[\begin{array}{cc}\mathrm{1}& \mathrm{0}\\ \mathrm{0}& -\frac{\mathrm{1}}{K\mathrm{\Delta }t}\\ \mathrm{0}& \mathrm{0}\end{array}\right],\end{array}$$

which can be evaluated numerically after convergence of the local iterations. The coefficients of this matrix represent the numerical values of the correspondent derivatives. They can be directly used in the evaluation of the tangent operator.

Next, we differentiate the expression for the volumetric-deviatoric decomposition of the Cauchy stresses (Eq. 3) with respect to the strain rate tensor, which yields:

$$\begin{array}{}\text{(B2)}& {\displaystyle \frac{\partial {\mathit{\sigma }}_{ij}}{\partial {\dot{\mathit{ϵ}}}_{kl}}}={\displaystyle \frac{\partial {\mathit{\tau }}_{ij}}{\partial {\dot{\mathit{\epsilon }}}_{mn}}}{\displaystyle \frac{\partial {\dot{\mathit{\epsilon }}}_{mn}}{\partial {\dot{\mathit{ϵ}}}_{kl}}}-{\mathit{\delta }}_{ij}{\displaystyle \frac{\partial p}{\partial {\dot{\mathit{\epsilon }}}_{\mathrm{II}}}}{\displaystyle \frac{\partial {\dot{\mathit{\epsilon }}}_{\mathrm{II}}}{\partial {\dot{\mathit{\epsilon }}}_{mn}}}{\displaystyle \frac{\partial {\dot{\mathit{\epsilon }}}_{mn}}{\partial {\dot{\mathit{ϵ}}}_{kl}}}.\end{array}$$

The derivative of the deviatoric strain rate (Eq. 4) with respect to the total strain rate gives the unit deviatoric projection tensor, which is given by:

$$\begin{array}{}\text{(B3)}& \begin{array}{rl}& {\displaystyle \frac{\partial {\dot{\mathit{\epsilon }}}_{ij}}{\partial {\dot{\mathit{ϵ}}}_{kl}}}=I_{ijkl}^D=I_{ijkl}-{\displaystyle \frac{\mathrm{1}}{\mathrm{3}}}{\mathit{\delta }}_{ij}{\mathit{\delta }}_{kl},\\ & I_{ijkl}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left({\mathit{\delta }}_{ik}{\mathit{\delta }}_{jl}+{\mathit{\delta }}_{il}{\mathit{\delta }}_{jk}\right).\end{array}\end{array}$$

By differentiating the norm of the effective deviatoric tensor with respect to the tensor itself, we obtain the scaled normalized deviatoric direction tensor given by the Eq. (49):

$$\begin{array}{}\text{(B4)}& {\displaystyle \frac{\partial {\dot{\mathit{\epsilon }}}_{\mathrm{II}}}{\partial {\dot{\mathit{\epsilon }}}_{kl}}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{n}_{kl}.\end{array}$$

The derivative of the deviatoric stress update expression (Eq. 48) reads:

$$\begin{array}{}\text{(B5)}& {\displaystyle \frac{\partial {\mathit{\tau }}_{ij}}{\partial {\dot{\mathit{\epsilon }}}_{kl}}}={\displaystyle \frac{\partial {\mathit{\tau }}_{\mathrm{II}}}{\partial {\dot{\mathit{\epsilon }}}_{\mathrm{II}}}}{\displaystyle \frac{\partial {\dot{\mathit{\epsilon }}}_{\mathrm{II}}}{\partial {\dot{\mathit{\epsilon }}}_{kl}}}{n}_{ij}+{\mathit{\tau }}_{\mathrm{II}}{\displaystyle \frac{\partial n_{ij}}{\partial {\dot{\mathit{\epsilon }}}_{kl}}}.\end{array}$$

Finally, differentiating the expression for the deviatoric direction tensor (Eq. 49) and using Eq. (B4) gives:

$$\begin{array}{}\text{(B6)}& {\displaystyle \frac{\partial n_{ij}}{\partial {\dot{\mathit{\epsilon }}}_{kl}}}={\displaystyle \frac{\mathrm{1}}{{\dot{\mathit{\epsilon }}}_{\mathrm{II}}}}{\displaystyle \frac{\partial {\dot{\mathit{\epsilon }}}_{ij}}{\partial {\dot{\mathit{\epsilon }}}_{kl}}}-{\displaystyle \frac{{\dot{\mathit{\epsilon }}}_{ij}}{{\dot{\mathit{\epsilon }}}_{\mathrm{II}}^{\mathrm{2}}}}{\displaystyle \frac{\partial {\dot{\mathit{\epsilon }}}_{\mathrm{II}}}{\partial {\dot{\mathit{\epsilon }}}_{kl}}}={\displaystyle \frac{\mathrm{1}}{{\dot{\mathit{\epsilon }}}_{\mathrm{II}}}}\left(I_{ijkl}-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{n}_{ij}{n}_{kl}\right).\end{array}$$

Substituting Eqs. (B3)–(B6) into Eq. (B2), making use of the matrix coefficients from Eq. (B1), and simplifying the resulting expression we get:

$$\begin{array}{}\text{(B7)}& {\displaystyle \frac{\partial {\mathit{\sigma }}_{ij}}{\partial {\dot{\mathit{ϵ}}}_{kl}}}=\mathrm{2}{\mathit{\eta }}_{\mathrm{eff}}{I}_{ijkl}^D+{\mathit{\beta }}_{\mathrm{1}}{n}_{ij}{n}_{kl}+{\mathit{\beta }}_{\mathrm{2}}{\mathit{\delta }}_{ij}{n}_{kl},\end{array}$$

where the effective stiffness constants are given by:

$$\begin{array}{}\text{(B8)}& {\mathit{\eta }}_{\mathrm{eff}}={\displaystyle \frac{{\mathit{\tau }}_{\mathrm{II}}}{\mathrm{2}{\dot{\mathit{\epsilon }}}_{\mathrm{II}}}},{\mathit{\beta }}_{\mathrm{1}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{A}_{\mathrm{11}}-{\mathit{\eta }}_{\mathrm{eff}},{\mathit{\beta }}_{\mathrm{2}}=-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{A}_{\mathrm{21}}.\end{array}$$

The derivative of the Cauchy stress with respect to the global pressure reads:

$$\begin{array}{}\text{(B9)}& {\displaystyle \frac{\partial {\mathit{\sigma }}_{ij}}{\partial p^{*}}}={\displaystyle \frac{\partial {\mathit{\tau }}_{ij}}{\partial p^{*}}}-{\mathit{\delta }}_{ij}{\displaystyle \frac{\partial p}{\partial p^{*}}}.\end{array}$$

Differentiating the deviatoric stress update expression (Eq. 48) and keeping in mind that the deviatoric direction tensor does not depend on pressure, yields:

$$\begin{array}{}\text{(B10)}& {\displaystyle \frac{\partial {\mathit{\tau }}_{ij}}{\partial p^{*}}}={\displaystyle \frac{\partial {\mathit{\tau }}_{\mathrm{II}}}{\partial p^{*}}}{n}_{ij}.\end{array}$$

After substituting Eq. (B10) into Eq. (B9) and using the coefficients from Eq. (B1) we obtain:

$$\begin{array}{}\text{(B11)}& {\displaystyle \frac{\partial {\mathit{\sigma }}_{ij}}{\partial p^{*}}}={\mathit{\beta }}_{\mathrm{3}}{n}_{ij}+{\mathit{\beta }}_{\mathrm{4}}{\mathit{\delta }}_{ij},\end{array}$$

where the effective stiffness constants are given by:

$$\begin{array}{}\text{(B12)}& {\mathit{\beta }}_{\mathrm{3}}=A_{\mathrm{12}},{\mathit{\beta }}_{\mathrm{4}}=-A_{\mathrm{22}}.\end{array}$$

For the visco-elastic case, A₁₁ becomes the only non-trivial matrix coefficient in the Eq. (B1). A quick inspection of Eq. (47) along with Eq. (B1) reveals that it can be directly evaluated as:

$$\begin{array}{}\text{(B13)}& A_{\mathrm{11}}={\left(A_{\mathrm{L}}+n{A}_{\mathrm{N}}{\mathit{\tau }}_{\mathrm{II}}^{n-\mathrm{1}}\right)}^{-\mathrm{1}}.\end{array}$$