Thibault Duretz
email: thibault.duretz@univ-rennes1.fr

René De Borst

Laetitia Le Pourhiet

On finite thickness of shear bands in frictional viscoplasticity, and implications for lithosphere dynamics

Keywords:

Numerical modeling of shear banding with

Introduction

Shear localization refers to the phenomenon of the concentration of strains in narrow zones when the applied load exceeds a certain threshold level. It occurs in virtually all materials [START_REF] Nadai | Plasticity[END_REF], takes place at all spatial and temporal scales within the lithosphere and manifests itself through phenomena which are widespread in rocks, e.g., faults, shear zones and shear bands. The localization of strains is observed in lithospheric domains where the stress levels are controlled mainly by either the temperature and the strain rate (viscous creep, ductile mode) or the pressure (frictional plasticity, brittle mode).

While the mechanisms that govern strain localization in the ductile mode are still being discussed [START_REF] Bercovici | A two-phase model for compaction and damage: 1. General theory[END_REF][START_REF] Précigout | Mantle weakening and strain localization: Implications for the long-term strength of the continental lithosphere[END_REF][START_REF] Thielmann | Shear heating induced lithospheric-scale localization: Does it result in subduction?[END_REF]Duretz and Schmalholz, 2015], those acting in the brittle mode, particularly in the realm of non-associated plasticity, are fairly well established [START_REF] Rudnicki | Conditions for the localization of deformation in pressure sensitive dilatant materials[END_REF][START_REF] Poliakov | Fractal plastic shear bands[END_REF].

For the purpose of geodynamic modeling, accounting for plasticity is mandatory to capture the self-consistent generation of tectonic plates [e.g. [START_REF] Tackley | Self-consistent generation of tectonic plates in time-dependent, threedimensional mantle convection simulations 2. strain weakening and asthenosphere[END_REF] and the development of fault zones [e.g. [START_REF] Poliakov | Initiation of salt diapirs with frictional overburdens: numerical experiments[END_REF][START_REF] Gerbault | Prediction of faulting from the theories of elasticity and plasticity: what are the limits?[END_REF]]. However, the physical processes within shear bands, which control their width, are often considered to be beyond the scope of current geodynamic models. Hence, for sake of simplicity, strain localization is often induced by a priori defined strain-softening functions [START_REF] Lavier | Self-consistent rolling-hinge model for the evolution of large-offset low-angle normal faults[END_REF][START_REF] Buck | A tale of two kinds of normal fault: The importance of strain weakening in fault development[END_REF][START_REF] Huismans | Asymmetric lithospheric extension: The role of frictional plastic strain softening inferred from numerical experiments[END_REF][START_REF] Buiter | The numerical sandbox: comparison of model results for a shortening and an extension experiment[END_REF][START_REF] Döhmann | Strain localization and weakening processes in viscously deforming rocks: Numerical modeling based on laboratory torsion experiments[END_REF], which are meant to take into account the role of complex thermohydro-chemico-mechanical interactions within faults in a phenomenological sense.

Since geodynamic models need to deal primarily with large deviatoric static strains, most models rely on the incompressible Stokes equations [e.g. [START_REF] Fullsack | An arbitrary lagrangian-eulerian formulation for creeping flows and its application in tectonic models[END_REF][START_REF] Tackley | Self-consistent generation of tectonic plates in time-dependent, threedimensional mantle convection simulations 2. strain weakening and asthenosphere[END_REF][START_REF] Gerya | Characteristics-based marker-in-cell method with conservative finite-differences schemes for modeling geological flows with strongly variable transport properties[END_REF]. Typically, visco-plasticity is used to describe the rheological behavior of geomaterials under these conditions, treating them as highly viscous fluids which can locally undergo plastic flow if a yield criterion is met [Willett, 1992]. Many geodynamic models neglect the role of elasticity and can therefore not capture the effects of elastic unloading, which can be pivotal in the proper description of progressive shearband development. Multi-dimensional stress states in geodynamics are typically captured using pressure-dependent yield functions like Drucker-Prager or Mohr-Coulomb [START_REF] Moresi | Incompressible viscous formulations for deformation and yielding of the lithosphere[END_REF][START_REF] Popov | SLIM3D: A tool for three-dimensional thermomechanical modeling of lithospheric deformation with elasto-visco-plastic rheology[END_REF].

For most rocks the angle of internal friction is relatively large in the steady state, around 30°or more [START_REF] Byerlee | Friction of rocks[END_REF]. On the other hand, the magnitude of the dilatancy angle is, which controls the amount of plastic volumetric change for a given amount of plastic shearing, is usually much lower, at most 10 -15°at the onset of shear banding [Vermeer and de Borst, 1984], tending to zero for progressively increasing deformations.

The latter observation ties in with the common use of the incompressible Stokes equations for these applications. The nearly incompressible nature of the plastic flow in combination with the strong frictional character of geomaterials renders the plasticity models non-associated.

Both the introduction of strain softening and the use of non-associated flow rules in pressure-dependent plasticity models yield mesh-dependent results and often exhibit an er-ratic and unsatisfactory convergence behavior of the equilibrium-searching iterative procedure [START_REF] Spiegelman | On the solvability of incompressible stokes with viscoplastic rheologies in geodynamics[END_REF]. The underlying reason is that the introduction of cohesion softening or of a non-associated flow rule has a mechanically destabilising influence. Indeed, both for cohesion softening [START_REF] Read | Strain softening of rock, soil and concrete -a review article[END_REF] and for non-associated flow [START_REF] Rudnicki | Conditions for the localization of deformation in pressure sensitive dilatant materials[END_REF] not only loss of mechanical stability can be induced, but worse, also loss of ellipticity. This mathematical condition is the basic cause of the frequently observed mesh dependence, which occurs in computations of materials which have constitutive laws that are equipped with this kind of behaviour. It has also been shown that under such conditions convergence of the incremental-iterative solution procedure deteriorates with increasing mesh refinement [START_REF] De Borst | Non-Linear Finite Element Analysis of Solids and Structures[END_REF].

For non-associated flow, loss of ellipticity can happen even if the material is still hardening and simulations have shown that global structural softening can then take place [START_REF] De Borst | Bifurcations in finite element models with a non-associated flow law[END_REF][START_REF] Pourhiet | Strain localization due to structural softening during pressure sensitive rate independent yielding[END_REF][START_REF] Sabet | Structural softening, mesh dependence, and regularisation in non-associated plastic flow[END_REF]. It occurs at more and more positive hardening rates when the difference increases between the angles of internal friction and dilatancy.

It is emphasized that the fundamental, mechanical-mathematical cause of these numerical problems is loss of ellipticity. Loss of ellipticity can cause the initial value problem to become ill-posed, which makes that solutions no longer continuously depend on the initial and boundary conditions. Numerical solutions then become meaningless, since they are fully dependent on the discretization, with respect to the fineness of the mesh, but also with respect to the direction of the grid lines [START_REF] Sluys | Discontinuous failure analysis for mode-i and mode-ii localization problems[END_REF][START_REF] Jirašek | Evaluation of directional mesh bias in concrete fracture simulations using continuum damage models[END_REF]. This holds for any discretization method, including meshless methods [START_REF] Pamin | Two gradient plasticity theories discretized with the element-free Galerkin method[END_REF], and also adaptive mesh refinement is severely biased [START_REF] Perić | On error estimates and adaptivity in elastoplastic solids: Applications to the numerical simulation of strain localization in classical and Cosserat continua[END_REF].

Unless a plasticity model which incorporates cohesion softening or non-associated flow is equipped with an internal length scale, it will lose ellipticity, and hence suffer from mesh sensitivity, at some loading stage. Yet, most plasticity models which have been used so far in long term tectonics do not incorporate an internal length scale. Such models, while ubiquitous in geophysics, geomechanics and engineering, are based on the assumption that the mechanical behavior in a point is representative for a small, but finite volume surrounding it. This assumption is often correct, but fails for highly localized deformations, like fault movement or shear bands. In the presence of strain weakening or non-associated flow local stress-strain relations have to be enriched to properly take into account the physical processes which occur at small length scales. A range of possibilities has been proposed to remedy this deficiency [START_REF] De Borst | Fundamental issues in finite element analysis of localisation of deformation[END_REF], including Cosserat plasticity [START_REF] Mühlhaus | The thickness of shear bands in granular materials[END_REF][START_REF] Stefanou | Cosserat approach to localization in geomaterials[END_REF], non-local plasticity [START_REF] Bažant | Nonlocal integral formulations of plasticity and damage: Survey of progress[END_REF] and gradient plasticity [START_REF] De Borst | Gradient-dependent plasticity: Formulation and algorithmic aspects[END_REF]. Unfortunately, all these possible solutions come with disadvantages. Obviously, they share the need to specify additional boundary conditions, which are often not easily understood from the physics. Other disadvantages are the introduction of additional degrees of freedom, as in Cosserat or gradient models, or they can result in fully populated, non-banded and non-symmetric stiffness matrices (non-local models).

For geodynamic applications the inclusion of a deformation-limiting viscosity, which has been tailored for modeling the deformation of crystalline solids [START_REF] Peirce | Material rate dependence and localized deformation in crystalline solids[END_REF][START_REF] Needleman | Material rate dependence and mesh sensitivity in localization problems[END_REF] represents an alternative to non-local rheological models. It is emphasized though, that not all visco-elasto-plastic rheologies solve the issue of mesh dependence, and that a pure series arrangement of the rheological elements (Maxwell-type approach, see Fig. 1b)) for instance, does not introduce a length scale and therefore does not remove the mesh-dependence issue . By contrast, a viscoplastic model which relies on the introduction of a rate-limiting viscosity in a parallel arrangement with a plastic slider (Kelvin-type viscoplasticity, Fig. 1c)) [START_REF] Perzyna | Fundamental problems in viscoplasticity[END_REF] does introduce a length scale and can provide mesh-independent numerical solutions [START_REF] Sluys | Wave propagataion and localization in a ratedependent cracked medium -model formulation and one-dimensional examples[END_REF]Wang et al., 1996;[START_REF] Dias Da Silva | A simple model for viscous regularization of elasto-plastic constitutive laws with softening[END_REF][START_REF] Niazi | Viscoplastic regularization of local damage models: Revisited[END_REF]. Although a Kelvin-type viscoplasticity rheology has been used in tectonic modeling studies before, [e.g. [START_REF] Hansen | A meshless formulation for geodynamic modeling[END_REF][START_REF] Regenauer-Lieb | The dynamics of multiscale, multiphysics faults: Part I -The long-term behaviour of the lithosphere[END_REF]Yin et al., 2018], the consequences have largely remained unexplored.

Herein, we will numerically study strain localization using non-associated Drucker-Prager plasticity, which captures the first-order behavior of the frictional lithosphere [START_REF] Lemiale | Shear banding analysis of plastic models formulated for incompressible viscous flows[END_REF][START_REF] Moresi | Incompressible viscous formulations for deformation and yielding of the lithosphere[END_REF][START_REF] Kaus | Factors that control the angle of shear bands in geodynamic numerical models of brittle deformation[END_REF]. We first illustrate the problem of mesh dependence using rate-independent plasticity. Then, we introduce a Kelvin-type rate-dependent viscoplastic formulation and demonstrate that the computed shear bands are mesh independent, even when strain softening is also introduced. We analyze the pressure, strain, and strain rate profiles across the shear bands as well as their evolution. Finally, we discuss the implications of using viscoplasticity for modeling in geodynamics.

In the remainder, we consider three rheological models. The first model is a standard, rate-independent elasto-plastic model (E-P), depicted in Fig. 1a) and, assuming small strains, it can be characterized by an additive decomposition of the strain rate into an elastic component and a plastic component:

= e + p , (1)
where the subscripts e and p denote elastic and plastic components, respectively. The deviatoric elastic strain e relates to the deviatoric stress τ τ τ as follows:

e = τ τ τ 2G , (2)
where G represents the shear modulus, which is kept constant for simplicity.

Plastic deformations arise when the yield function

F = J 2 -C cos(φ) -P sin(φ), (3)
attains a zero value, with C and φ the cohesion and the angle of internal friction, respectively.

J 2 = 1 2 τ 2 xx + τ 2 yy + τ 2 zz + τ 2
xy is the second invariant of the deviatoric stresses τ xx , τ yy , τ zz , τ xy and P = -1 3 σ xx + σ yy + σ zz is the mean stress, defined as negative in tension. When cohesion hardening or softening is incorporated, the hardening/softening modulus takes the form:

h = dC d p , ( 4
)
where p is the accumulated equivalent plastic strain according to the strain-hardening hypothesis and is formulated as:

p = ∫ 2 3 ( p ) T p dt. ( 5
)
During continued plastic flow, F = 0 and the deviatoric strain rates are assumed to be derivable from a plastic potential function Q:

p = λ ∂Q ∂σ σ σ , ( 6
)
where λ is a plastic multiplier and Q is assumed to have a form similar to that of the yield function F:

Q = J 2 -P sin(ψ), (7)
with ψ ≤ φ the dilation angle.

Next, we consider an visco-elasto-plastic model classically used in geodynamics, Fig. 1b), with a viscous damper added in series to the previous rheology (Maxwell-type approach, V-E-P):

= v + e + p , (8)
the superscript v denoting a viscous (strain) component. The viscous deviatoric strain rate v is assumed to be linearly related to the deviatoric stress tensor:

v = τ τ τ 2η (9)
with η the (constant) dynamic shear viscosity.

Thirdly, we consider a model where the plastic element of Fig. 1a) is substituted by a viscoplastic element, Fig. 1c), which can be considered as a Kelvin-type arrangement (E-VP):

= e + vp . (10)
During viscoplastic flow, the yield function is now defined as [START_REF] Heeres | A comparison between the Perzyna viscoplastic model and the Consistency viscoplastic model[END_REF]:

F = J 2 -C cos(φ) -P sin(φ) -η vp λ (11)
where η vp is the viscosity of the damper. The rate-independent limit is recovered by letting η vp → 0. Expression (11) makes the yield function rate-dependent, so this model belongs to the class of consistency viscoplastic models. It has been shown [Wang et al., 1997;[START_REF] Heeres | A comparison between the Perzyna viscoplastic model and the Consistency viscoplastic model[END_REF] that this class of viscoplastic models has advantages over overstresstype viscoplastic models, e.g. those of the Perzyna-type [START_REF] Perzyna | Fundamental problems in viscoplasticity[END_REF], including an improved convergence behaviour and a more straightforward implementation.

Numerical Implementation

The expression of the visco-elastic tangent matrix D ve is obtained by integrating the Maxwell rheological chain, Eq. ( 8) under the assumption of no plastic flow. We introduce the quantities

G ve = 1 G + 1 η -1 and ξ = G ve G (12)
to obtain the following update rule for the total stress tensor:

σ σ σ t+1 = -P t + ξτ t + D ve ∆ ve t+1 (13)
where the ∆-operator represents a finite increment, and the visco-elastic tangent operator reads:

D ve =               K + 4 3 G ve K -2 3 G ve K -2 3 G ve 0 K -2 3 G ve K + 4 3 G ve K -2 3 G ve 0 K -2 3 G ve K -2 3 G ve K + 4 3 G ve 0 0 0 0 G ve               , (14)
with K the elastic bulk modulus.

If (visco)plastic flow has occurred, the incremental plastic multiplier, ∆λ must be computed from Eq. ( 11) with λ = ∆λ ∆t and F = 0, see also [START_REF] Heeres | A comparison between the Perzyna viscoplastic model and the Consistency viscoplastic model[END_REF]. Using a

Taylor's expansion for the yield function [START_REF] De Borst | Studies in anisotropic plasticity with reference to the Hill criterion[END_REF][START_REF] Duretz | The benefits of using a consistent tangent operator for viscoelastoplastic computations in geodynamics[END_REF],

or by considering that the corrected stress state lie onto the yield surface [START_REF] De Souza Neto | Computational Methods for Plasticity: Theory and Applications[END_REF], a closed-form expression for ∆λ can be derived for a Drucker-Prager yield function:

∆λ = F(σ σ σ trial ) G ve + K sin(φ) sin(ψ) + η vp ∆t + H . ( 15
)
where σ σ σ trial is the trial stress, which has been computed assuming no (visco)plastic flow, and

H = h cos φ 2 3 ∂Q ∂σ σ σ T ∂Q ∂σ σ σ . ( 16
)
Defining σ σ σ t as the stress state at the beginning of the loading step, the new stress state can be computed by adding the viscoelastic stiffness times the difference of the total and the (visco)plastic strain increments to σ σ σ t :

σ σ σ t+1 = σ σ σ t + D ve (∆ -∆ vp ) = σ σ σ t + D ve ∆ -∆λ ∂Q ∂σ σ σ (17)
which exactly satisfies the yield condition.

When a Newton-Raphson iterative procedure is used to achieve global equilibrium, the above expression must be linearized. This leads to the so-called consistent tangent operator for visco-elastic-viscoplastic solids

D vep ≡ ∂σ σ σ ∂ = E -1 D ve - E -1 D ve ∂Q ∂σ σ σ ∂F ∂σ σ σ T E -1 D ve H + η vp ∆t + ∂F ∂σ σ σ T E -1 D ve ∂Q ∂σ σ σ (18)
with ∆t the time step and

E = I + ∆λ ∂ 2 Q ∂σ σ σ 2 . ( 19
)
Detailed derivations of the consistent tangent operator and incremental plastic multiplier are provided in the appendix.

Model configuration

The results shown in this study have been obtained using a simple initial model configuration. The two-dimensional model consists of a 1.0 × 0.7 domain subjected to a kinematic boundary condition which induces a pure shear state. Displacements increments (∆u BC i = x i ∆ BC ) are imposed on the south and east sides of the domain. The west and north sides of the domain are slip-free boundaries. A circular inclusion with radius 5 × 10 -2 is located at the origin of the domain. This imperfection is characterized by a lower shear modulus, which causes the stress perturbation ultimately leading to strain localization. All initial stress and strain components are set equal to zero. We have used the same number of nodes (n i ) in both spatial dimensions (n x = n y ). The shear modulus G has been set equal to 1 in the matrix and equal to 0.25 in the perturbation, while the bulk modulus K has been set equal to 2. The applied strain increment ∆ BG = 5.0 × 10 -6 . For the viscoplastic model, the viscosity is set to 2.5 × 10 2 and the timestep ∆t = 10 4 , which yields a background strain rate BG = ∆ BG ∆t = 5.0×10 -9 . More information about the model parameters is given in Table 1.

Modeling results with an elasto-plastic rheology (E-P model)

The first series of computations have been carried using a rate-independent elastoplastic rheology, an angle of internal friction φ = 30 o and a dilatancy angle ψ = 10 o .

Three different resolutions were employed, with n i = [51, 101, 201] nodes. Fig. 2 shows the second invariant of the accumulated strain for all three resolutions, for the same amount of applied background strain. A single shear band develops starting from the imperfection.

The shear band is oriented at 35 o from the direction of the principal compressive stress, which is in line with the Arthur formula (45 o -1/4(φ + ψ)) [START_REF] Arthur | Plastic deformation and failure in granular media[END_REF][START_REF] Kaus | Factors that control the angle of shear bands in geodynamic numerical models of brittle deformation[END_REF], which has been experimentally observed for shear banding in sands and can be derived from bifurcation analysis using the Mohr-Coulomb criterion. For the Drucker-Prager yield criterion the out-of-plane stress, however, plays a role, [START_REF] Rudnicki | Conditions for the localization of deformation in pressure sensitive dilatant materials[END_REF], but this apparently affects the numerical results only marginally. The results are clearly mesh dependent as the localized strain is distributed over a thickness of a single cell.

Another representation was made by plotting the profiles of the second invariant of the accumulated strain invariant and of the pressure across the shear bands, see Fig. 2d). The profiles reveal a divergence of the solutions with increasing resolution. Since the displacement jump across the shear band is constant, the magnitude of strain in the shear band continues to increase for finer resolutions. For the finest discretization we observe the development of a Dirac-like strain distribution. This is further evidence for the illposed character of the boundary value problem, since the Dirac-like strain profile points at a discontinuity in the displacements, which can only occur if the governing equations have locally changed character from elliptic to hyperbolic.

Results for an elasto-viscoplastic rheology (E-VP model)

The second series of models were carried out using a elasto-viscoplastic rheology.

The numerical simulations were achieved on progressively refined meshes consisting of [51,101,201,401,801] nodes. Fig. 3 depicts the spatial distribution of the accumulated strain at a fixed amount of shortening (3.0 × 10 -4 ). In contrast with the elasto-plastic models of Fig. 2 shear bands of a finite width now arise. We note that the shear band is still oriented at 35 o from the direction of σ 1 . Detailed probing of different variables, in particular the pressure, the second invariant of the accumulated strain and strain rate, results in a clear convergence upon mesh refinement, see Fig. 4. The strain and strain rate profiles across the shear bands have a quasi-Gaussian shape. The peak strain and the peak strain rate are at the center of the shear band and reach values of 1.8×10 -3 and 5.3×10 -9 , respectively. The pressure is lower inside the shear band and reaches a minimum value of 0.9 × 10 -4 . For a low resolution (51 2 nodes), the peak strain reaches about half the magnitude of that obtained with the higher resolutions and the strain is localized over a wider zone.

n i =
The evolution of the strain, the strain rate and the pressure across the shear band is shown in Fig. 5a) for a given fixed resolution. From a background strain of 1.0 × 10 -4 , the strain locally increases inside a well-defined region. The amplitude of a Gaussian-like strain profile exhibits a growth of 1.8 × 10 -3 over an increment of background strain of 2.0 × 10 -4 , see Fig. 5a). The evolution of the second invariant of the strain rate shares these characteristics. The amplitude of the Gaussian-like profile of the strain rate reaches a peak value of 6.2 × 10 -9 shortly after the onset of shear localization (background strain 2.0 × 10 -4 ). This corresponds to a magnitude of the strain rate which is approximately twelve times larger than the applied background strain rate. The profiles of the pressure are characterized by a progressive decrease towards the center of the shear zone. After an applied background strain of 3.0 × 10 -4 , the pressure in the shear band drops to about two-thirds the background pressure. The shape of the pressure profile differs from that of the accumulated strain rate profiles. With ongoing strain, it progressively deviates from a Gaussian shape and morphs into a square shape. At the final stage, positive pressures develop adjacent to the shear zone, which further increases the magnitude of the pressure gradient across the shear band.

Characteristic shear band width

The above results indicate that the viscoplastic rheology introduces a length scale into the boundary value problem. In the following, we define the characteristic shearband thickness (D vp ) as the bandwidth of the accumulated strain profiles, which can be extracted from the two-dimensional modeling. Since the strain rate profiles sampled normal to the shear band exhibit a Gaussian-like shape, they can approximated as:

= max exp - z 2 D vp 2 , ( 20
)
where z is coordinate orthogonal to the shear band and max is the maximum value of second strain rate invariant along the profile. The value of D vp is obtained by an optimum fit of the Gaussian equation to the profiles which have been extracted from the twodimensional models.

The transient evolution D vp for different mesh resolutions is given in Fig. 6a). For all resolutions, the initial value of D vp is equal to the radius of the initial perturbation (r = 0.05). The values of D vp progressively decrease with an increasing applied strain. They all reach an asymptotic value at an applied strain of approximately 2.5 × 10 -4 . We have plotted the asymptotic D vp values as a function of the grid spacing h -1 in Fig. 6b). The shear band width asymptotically approaches a value of 0.01 with an increasing resolution.

Using a similar rheological model, Wang et al. [1996] have quantified the impact of the dimension of the initial imperfection on the shear band width. The dimensions of the initial perturbation are, together with the material parameters, key to the occurrence and further evolution of shear localization. It plays a fundamental role at the onset of shear localization. However, when the shear zone reaches a steady-state situation, i.e. when the width D vp has stabilized, the shear band widths are virtually independent of the size of the imperfection, as shown in Fig. 7). This is similar to results obtained with thermomechanical models of strain localization in temperature activated rate-dependent materials [START_REF] Lemonds | Finite element analyses of shear localization in rate and temperature dependent solids[END_REF]] using a power law viscous rheology [START_REF] Duretz | Physics-controlled thickness of shear zones caused by viscous heating: Implications for crustal shear localization[END_REF] and a power law visco-elastic rheology [Duretz et al., 2015]. In a steady state, the characteristic shear band thickness is essentially independent of the size of the imperfection.

Modeling strain localization at the crustal scale

In order to investigate the applicability and consequences of viscoplasticity for modeling lithospheric deformations, we have carried out simulations on kilometer-scale dimensions, using typical material parameters for rocks (Model Crust 1 and Model Crust 2, see Table 1). Model Crust 1 was designed to study the initiation and propagation of a single shear band which originates from a well-defined material imperfection using cohesion softening. The configuration shares similarities with that of Model E-VP, but has a material imperfection with radius of 500 m. The shear modulus within the perturbation is equal to 25% of that in the matrix (Table 1). The boundary velocity was set to

V BC = 2 × 10 -9 m.s -1
. The dimensions of the domain are 10 × 6.85 km, which is discretized using 400 2 cells leading to a resolution of 25 m. For the reference test, the Kelvin element viscosity was set to 10 17 Pa • s, the initial cohesion was set to C = 1.75 × 10 7 Pa.

Softening was prescribed by setting a negative value of hardening modulus, h = -7.0 × 10 7

Pa and allowing a reduction of cohesion by a factor 2. The shear band develops from the south-west corner towards the north-east corner, see Fig. 8a). As for the previous cases, strain localization is progressive and the shear band width narrows down with increasing time or strain, Fig. 8b), and progressively reaches a width of 90 m. In order to test the sensitivity of the model, we have tested different parameter combinations, all satisfying a constant value the product η vp V BC C . The resulting models all predict a final shear band thickness of about 90 m. The time needed for strain localization is proportional to the Kelvin element viscosity. Defining the characteristic time t c ∝ η vp G and the characteristic

length D c ∝ η vp V BC C
allows to collapse the shear band thickness evolution onto a single master curve (Fig. 8c).

When considering crustal scale strain localization, incompressible plastic deformation is generally invoked. However, with rate-independent plastic models, such a limit poses serious numerical issues. The latter occur for large differences between the friction and dilation angles (φψ > 20 o ), and are caused by the loss of ellipticity [START_REF] Sabet | Structural softening, mesh dependence, and regularisation in non-associated plastic flow[END_REF]. In practice, models often diverge (e.g. [START_REF] Spiegelman | On the solvability of incompressible stokes with viscoplastic rheologies in geodynamics[END_REF]) even when using consistent linearizations (e.g. [START_REF] Duretz | The benefits of using a consistent tangent operator for viscoelastoplastic computations in geodynamics[END_REF]). With Kelvin-type viscoplasticity models, attainment of convergence appears to be much less problematic because of the weak regularisation of the ill-posed problem. Fully converged results were obtained for ψ = 0 • and φ = 30 • as depicted on Fig. 9a. Interestingly, strain localisation was obtained for values of the dilation angle up to 10 • , but localisation did not occur for larger values of ψ. During localisation, the shear zone narrows down at a faster rate when the dilation angle is small (Fig. 9b). However, the width of the post-localisation shear zone is not affected by a variations of ψ (Fig. 9b). By contrast, stress-strain curves notably differ for the different values of φ. In the incompressible limit, the effective stress ( σ

= 1 V ∫ √ J II dV)
reaches a peak value (15 MPa) and then decreases to a saturation value (≈ 10 MPa). For larger angles, the effective stress keeps increasing despite the occurrence of strain localization.

In Model Crust 2 strain localization is seeded by setting an initial random perturbation on the cohesion field. The confinement pressure was set equal to 50 MPa and cohesion softening is again applied. Numerous intersecting shear bands of different lengths develop, see Fig. 10). Due to the complex internal kinematics, shear bands exhibit different lengths and widths. This is in contrast with previously presented models in which single shear bands were arising from from single perturbations (Fig. 8a). We have run simulations for various resolutions [101,201,401,801]) up to a bulk strain of 0.5 × 10 -2 .

(n i =
The results show that, despite the use of non-associated plasticity and cohesion softening, a reasonable convergence upon mesh refinement was achieved. While there is a difference in the strain fields between the low resolution models, e.g. between 100 2 and 200 2 cells, the differences are much less pronounced when comparing simulations for higher resolutions, e.g. between 400 2 and 800 2 cells. Most importantly, global equilibrium iterations converged quadratically up to machine precision without any need to reduce the applied strain increment. This is in complete contrast with rate-independent elasto-plastic models, with which it was impossible to reach such a high accuracy (results not shown here).

Discussion

A potential regularization for mesh-dependent strain softening plasticity

In tectonic modeling, strain localization in the frictional domain is generally modeled using a strain-softening parametrization, e.g. [START_REF] Huismans | Asymmetric lithospheric extension: The role of frictional plastic strain softening inferred from numerical experiments[END_REF]. A progressive decrease of the magnitude of plastic parameters (cohesion, friction angle) is imposed as a function of the accumulated plastic strain. Such a parametrization is supposed to mimick the effects of complex physico-chemical processes (e.g. metamorphic reactions, fluid pressure variations) that are not taken into account in the model description, which may cause post-localization stress drops measured in experiments on rock deformation.

However, this approach is known to suffer from mesh dependence, cf. Fig. 2: the results strongly depend on the numerical resolution. Among the numerous possibilities that provide regularization of mesh dependence in plasticity, viscoplasticity is a simple and efficient solution. Herein we have shown that using a viscoplastic rheology together with a non-associated frictional plasticity model and cohesion softening provides results which converge upon mesh refinement. Morerover, we have found that viscoplasticity also facilitates resolving strain localization in the incompressible plastic limit (ψ = 0), which is very relevant for lithospheric conditions.

Benefits for the convergence behavior of the global non-linear solver

Use of a viscoplastic model as regularization method can also dramatically improve global convergence and the computability and solvability of shear banding. With a viscoplastic model, strains are not concentrated in one cell, but are distributed over a finite width. Therefore, increasing the resolution does not lead to an increase of the strain locally, which can cause serious numerical issues such as local snap-backs in the return mapping and the occurrence of multiple, non-physical equilibrium states which cause divergence of the global equilibrium-finding iterative procedure. An illustration is that for an elasto-plastic rheology no results could be obtained for meshes with more than 201 2 nodes, even when reducing the strain increment. When the strain increment was kept constant, it was not possible to reach convergence for resolutions with more than 101 2 nodes (Table 2). Using an elasto-viscoplastic rheology, however, converged results were obtained for fine resolutions, up to 801 2 nodes (Figs 3 and 10). Moreover, both the maximum and average number of iterations required to achieve global equilibrium are almost insensitive to the numerical resolution, as expected (Table 2).

Implications of viscoplasticity

The viscoplastic model introduces a rate dependence in the yield function and requires an additional model parameter: the viscosity of the viscoplastic (Kelvin) element.

Here we used the viscosity as a numerical parameter rather than as an experimentally measured quantity with a clear physical meaning. Shear viscosities estimated from labora-tory experiments of rock deformation exhibit an Arrhenius temperature dependence, which results in a exponential growth of viscosity with decreasing temperature. Using such values as viscosity of the Kelvin element would lead to an unrealistic overshoot of the yield function in the frictional (low temperature) domains of the lithosphere and would preclude strain localization. We believe that the viscosity of the Kelvin element should be chosen such that shear bands can be resolved numerically. We expect that such an approach will allow for the converged resolution of shear bands in geodynamic models where the spatial resolution has so far at most reached the 100 meter scale. Another approach to select the Kelvin element viscosity is to reason in terms of the overstress (∆σ = η vp λ). By assuming that the rate of the plastic multiplier is proportional to the background strain rate ( λ ∝ LV BC ), it is possible to define the viscosity that will approximately generate the predefined overstress (η vp ∝ ∆σ LV BC ). For example, using the parameters of Model Crust 1 (Fig. 8) and assuming an overstress of 20 kPa, one obtains a viscosity for the Kelvin element of 10 17 Pa•s.

In general, the width of shear bands which arise in the frictional domains of the lithosphere are highly variable and can range from discrete fault planes to finite thickness fault zones (gouges), which involve complex processes that are beyond the scope of this study (shear heating, fluid pressure variations, grain crushing, mineral reactions). A detailed study of these processes may provide a physics-based regularization for the width of frictional shear bands, which will likely be smaller than the current resolution power of geodynamic models.

Differences with the standard visco-elasto-plastic model

In geodynamics, visco-elasto-plasticity is generally implemented via a rheological model which couples a viscous damper, an elastic spring and a plastic slider in series, which can be dubbed a Maxwell V-E-P model, e.g. [START_REF] Lemiale | Shear banding analysis of plastic models formulated for incompressible viscous flows[END_REF]; [START_REF] Gerya | Robust characteristics method for modelling multiphase visco-elasto-plastic thermo-mechanical problems[END_REF]; Kaus [2010]. With such a model, rate dependence is included in the visco-elastic trial stress, but not in the plastic strain component. Shear localization obtained with such models has the same characteristics as that obtained with a rate-independent elasto-plastic rheology. Shear bands localize on a single band of cells or elements, thus causing numerical simulations to be mesh sensitive, see Figs. 11a,b). This is in contrast with the model discussed above, which incorporates rate dependence in the plastic element and allows shear localization to spread over several cells or elements, see Fig. 11c).

Simulations with an elasto-plastic model may also require extremely fine load increment to reach global equilibrium, as shown in Fig. 11d). This equally applies to models with a viscoelastic-plastic (V-E-P) rheology. Elasto-viscoplastic (E-VP) models overcome this issue (Fig. 11d) and can accurately compute the load-bearing capacity of an E-P material, see Fig. 11e).

For geodynamics modeling purposes, we suggest that E-VP and V-E-P models could be combined into an E-V-VP model, shown in Fig. 1d). One one hand, Maxwell viscoelasticity is necessary to capture both the short timescale [START_REF] Deng | Viscoelastic flow in the lower crust after the 1992 landers, california, earthquake[END_REF][START_REF] Heimpel | Earthquake scaling: the effect of a viscoelastic asthenosphere[END_REF]Wang, 2007] and the long timescale, e.g. [START_REF] Farrington | The role of viscoelasticity in subducting plates[END_REF]; Schmalholz et al.

[2015]; [START_REF] Olive | The role of elasticity in simulating long-term tectonic extension[END_REF], which are essential features of lithospheric deformations. On the other hand, Kelvin-type viscoplastic models remedy known issues in modeling strain localization in the lithosphere. A combined E-V-VP model would be suitable to capture the visco-elastic behavior of rocks, but would also enable to obtain mesh-independent and globally convergent solutions of plastic shear banding.

Conclusions

We have investigated the role of elasto-viscoplasticity with a damper in parallel to a plastic slider (Kelvin-type rheology) on the development of shear bands in the frictional regions of the lithosphere. While the rate-independent frictional plasticity models, which have been used classically, suffer from mesh sensitivity, models using this viscoplastic rheology converge upon mesh refinement. The strain, the strain rate and the pressure inside the shear bands reach finite values upon a decrease of the grid spacing. A characteristic length scale is introduced due to the rate dependence of the viscoplastic model.

Our results indicate that shear bands which arise from pressure-dependent viscoplasticity maintain their orientation, but are now also equipped with a characteristic band width.

Even the combination of a Kelvin-type viscoplasticity with strain softening on the cohesion gives a mesh-convergent behavior. The approach is thus a viable way to regularize strain localization in geodynamic models. Most importantly, the introduction of an internal length scale due to viscoplasticity maintains well-posedness of the boundary value problem also during shear banding, and therefore markedly improves the convergence of equilibrium iterations, which is a recurrent issue in geodynamic simulations.

11 Figures Results were computed on a mesh consisting of 401 2 nodes.

A: The incremental plastic multiplier

For (visco)plastic deformations to occur, the stress must lie on the yield surface and continue to be on the yield surface for an infinitesimal time increment when adopting the consistency model of viscoplasticity [Wang et al., 1997;[START_REF] Heeres | A comparison between the Perzyna viscoplastic model and the Consistency viscoplastic model[END_REF]. This implies that for the yield function at the end of the plastic increment, we must have F t+1 = 0.

Following [START_REF] De Souza Neto | Computational Methods for Plasticity: Theory and Applications[END_REF], the corrected total stress is given by:

σ σ σ t+1 = 1 - G ve ∆λ J trial I I σ σ σ trial (A.1)
Using the identity σ σ σ t+1 √ J t+1

I I = σ σ σ trial √ J trial I I
, the corrected second stress invariant is:

J t+1 I I = J trial I I -G ve ∆λ. (A.
2)

The corrected pressure reads:

P t+1 = P trial + K sin(ψ)∆λ. (A.3)
and the updated cohesion can be expressed as:

C t+1 = C trial + H∆λ (A.4)
Approximating the rate of the plastic multiplier as λ = ∆λ ∆t , the yield function at the end of the plastic increment can be written explicitly as:

Figure 1 .

 1 Figure 1. Investigated rheological models for deformation of the lithosphere: a) Maxwell elasto-plastic model. b) Maxwell visco-elasto-plastic model. c) Elastic model coupled to a Kelvin viscoplastic element. d) Maxwell visco-elastic model coupled to a Kelvin viscoplastic element.

Figure 2 .

 2 Figure 2. Spatial distribution of accumulated strain ( I I ) calculated with an elasto-plastic rheology for three different mesh resolutions (51 2 , 101 2 , and 201 2 nodes). Results are depicted after a bulk strain of ≈ 7.7 × 10 -5 . The white lines indicate the location of solution profiles reported in d. d) Profiles of accumulated strain (a) probed across elasto-plastic shear bands.

Figure 3 .

 3 Figure 3. Spatial distribution of accumulated strain ( I I ) calculated with an elasto-viscoplastic rheology for five different mesh resolutions (51 2 , 101 2 , 201 2 , 401 2 and 801 2 nodes). Results are depicted after a bulk strain of 3.0 × 10 -4 . The white lines indicate the location of solution profiles reported in Fig. 4.

Figure 4 .

 4 Figure 4. Profiles of accumulated strain (a), effective strain rate (b) and pressure (c) probed across elastoviscoplastic shear bands. The results were obtained on five different mesh resolutions (51 2 , 101 2 , 201 2 , 401 2and 801 2 nodes). The solution profiles were sampled along the white lines visible on Fig.3

Figure 5 .

 5 Figure 5. Temporal evolution of accumulated strain, effective strain rate, and pressure across the shear zone. The different profiles correspond to five values of background strain (from 1.0 × 10 -4 to 3.0 × 10 -4 ).

Figure 6 .

 6 Figure 6. a) Evolution of the characteristic shear zone thickness (D vp ) versus accumulated background strain ( BG ). The values were extracted from the runs with five different mesh resolutions (Fig. 3). b) Characteristic shear band thickness (D vp ) versus grid spacing h. The results were extracted from the runs depicted on Fig. 3

Figure 7 .

 7 Figure 7. Influence of the initial size of the imperfection on the shear bands. Models were run with five different initial seed radius and for a fixed resolution of 401 2 nodes. The radius of the imperfection in the reference model is r (panel c). The inclusion has a reduced size in panels a) and b) and an increased size in panels d) and e). No material softening was applied.

Figure 8 .

 8 Figure 8. Strain localization at the kilometer scale arising from a single material perturbation. Panel a) depicts the strain field after a bulk strain of 5 × 10 -3 . Panel b) shows measured shear band widths for different parameter combinations. Panel c) depicts the master curve obtained when defining the characteristic time scale directly proportional to the Kelvin element viscosity.

Figure 9 .

 9 Figure 9. Effect of plastic dilatancy on shear band development. Panel a) depicts the strain field after a bulk strain of 1 × 10 -2 in the incompressible limit (ψ = 0 • ). Panel b) shows measured shear band widths for different values of dilatancy angle. Panel c) depicts stress-strain curves for the different values of ψ.

Figure 10 .

 10 Figure 10. Shear banding at kilometer scale arising from an initial random perturbation of the initial cohesion field for various mesh resolutions (from 100 2 to 800 2 cells). A Kelvin E-VP rheological model was applied. The confining pressure was set to 50 MPa and the Kelvin element viscosity was set to η VP = 3 × 10 18 Pa.s.

Figure 11 .

 11 Figure 11. Spatial distribution of accumulated strain ( I I ) after a bulk strain of 3.0×10 -4 . a) depicts results obtained with an elasto-plastic rheology (E-P). b) corresponds to a visco-elasto-plastic rheology (V-E-P) using a serial viscosity of 2.5 × 10 5 . and c) show results obtained with an elasto-viscoplastic rheology (E-VP) using a Kelvin viscosity of 2.5 × 10 2 . Panel d) shows the variations of the strain increment needed for achieving successful non-linear solutions. Panel d) shows the evolution of effective stress for the 3 different rheological models.

F

  t+1 = J trial I I -G ve ∆λ -cos(φ) C trial + H∆λ -sin(φ) P trial + K sin(ψ)∆λ -η vp ∆λ ∆t . (A.5)Solving for F t+1 = 0 then yields the following expression for the incremental plastic multiplier:∆λ = F trial G ve + K sin(φ) sin(ψ) + η vp ∆t + H , (A.6)where F trial = J trial I I -sin(φ)P trialcos(φ)C trial .B: The visco-elastic-plastic consistent tangent operatorDuring visco-elasto-(visco)plastic straining, the stress update follows:σ σ σ t+1 = -P t i + ξτ τ τ t + D ve ∆ -∆λD ve ∂Q ∂σ σ σ (B.1)A small variation δ of the updated stress σ σ σ t+1 is given by:δσ σ σ = D ve δ -δλD ve ∂Q ∂σ σ σ -∆λD ve ∂ 2 Q ∂σ σ σ 2 δσ σ σ. (B.2)and can be recast as:δσ σ σ = E -1 D ve δ -E -1 D ve ∂Q ∂σ F = F(σ σ σ, λ, λ), can be elaborated as: approximation δ λ = δλ ∆t , and invoking condition (B.6), provides an expression for the variation of the plastic multiplier: ∂F ∂λ and η vp ≡ ∂F ∂ λ . This expression is now substituted into Eq. (B.3), which leads to:δσ σ σ = E -1 D ve -

Table 1 .

 1 List of parameters relative to the different tests presented in this study. *: In Models 1 and 3, the strain increment may vary with time to allow for global convergence (see Fig.11d). **: in Model Crust 2, the initial cohesion is randomly perturbed.

	Parameter Model 1 E-P Model 2 E-VP Model 3 V-E-P Model 4 E-VP soft. Model Crust 1 Model Crust 2
	L x [m]	1	1	1	1	1.41 × 10 4	2 × 10 4
	L y [m]	0.7	0.7	0.7	0.7	0.685 × 10 4	10 4
	∆ BC [-]	5 × 10 -6 *	5 × 10 -6	5 × 10 -6 *	5 × 10 -6	2 × 10 -5	5 × 10 -5
	∆t [s]	-	10 4	10 4	10 4	10 8	10 10
	P c [Pa]	0	0	0	0	0	50 × 10 6
	C [Pa]	1.75 × 10 4	1.75 × 10 4	1.75 × 10 4	1.75 × 10 4	1.75 × 10 7	10 7 * *
	h [Pa]	0	0	0	-10 -2	-7 × 10 7	-1.5 × 10 7
	K [Pa]	2	2	2	2	2 × 10 10	2 × 10 10
	G [Pa]	1	1	1	1	10 10	10 10
	η [Pa.s]	-	-	2.5 × 10 5	-	-	-
	η VP [Pa.s]	-	2.5 × 10 2	-	2.5 × 10 5	1 × 10 17	3 × 10 18

Table 2 .

 2 Number of Newton-Raphson iterations required to reach global equilibrium for different resolutions using either the E-P or the E-VP model. The relative tolerance was set to 10 -11 and was measured using the L 2 norm. Runs which did not result in a converged state are denoted by a dash.

	E-P models	51 2 101 2 201 2 401 2
	Mean # its.	3.88 6.78	-	-
	Max. # its.	11	26	-	-
	E-VP models 51 2 101 2 201 2 401 2
	Mean # its.	3.88 4.17 4.41 4.89
	Max. # its.	11	12	12	12

c) Elasticity with Kelvin viscoplasticity (E-VP)

Confidential manuscript submitted to Geochemistry, Geophysics, Geosystems

Acknowledgments

The research reported in this article has been partially supported by the European Research Council under grant 664734 "PoroFrac". The data presented in this study are the result of numerical simulations. Results can be reproduced using the open source routines M2Di_EP (https://bitbucket.org/lraess/m2di/src/e4bde7404e533ffacac0746872c5b0d24593517b/ M2Di_EP/?at=master).

Confidential manuscript submitted to Geochemistry, Geophysics, Geosystems Confidential manuscript submitted to Geochemistry, Geophysics, Geosystems Vermeer, P. A., and R. de Borst (1984), Non-associated plasticity for soils, concrete and rock, Heron, 29(3), 3-64. Wang, K. (2007), Elastic and viscoelastic models of crustal deformation in subduction earthquake cycles, in The Seismogenic Zone of Subduction Thrust Faults, edited by T. Dixon and J. C. Moore, pp. 540-575, Columbia University Press, New York. Wang, W. M., L. J. Sluys, and R. de Borst (1996), Interaction between material length scale and imperfection size for localisation phenomena in viscoplastic media, European Journal of Mechanics. A, Solids, 15, 447-464.

Wang, W. M., L. J. Sluys, and R. de Borst (1997), Viscoplasticity for instabilities due to strain softening and strain-rate softening, International Journal for Numerical Methods in Engineering, 40(20), 3839-3864. Willett, S. D. (1992), Dynamic and kinematic growth and change of a Coulomb wedge, in Thrust Tectonics, edited by K. R. McClay, pp. 19-31, Springer, Dordrecht. Yin, A., Z. Xie, and L. Meng (2018), A viscoplastic shear-zone model for deep (15-50 km) slow-slip events at plate convergent margins, Earth and Planetary Science Letters,491,[81][82][83][84][85][86][87][88][89][90][91][92][93][94]