Multiphysics Modeling of a Brittle‐Ductile Lithosphere: 1. Explicit Visco‐Elasto‐Plastic Formulation and Its Numerical Implementation

The long‐term strength of the lithosphere is controled by two different modes of deformation: a brittle‐like, effective pressure‐sensitive behavior at shallow crustal depth, which gradually transits to a thermally activated ductile flow rheology with increasing depth. All applications dealing with long‐term tectonics therefore share the necessity to describe in a consistent way the multiphysics coupling among the different deformation mechanisms controlling the bulk behavior of the lithosphere. We describe an efficient numerical implementation of a consistent visco‐elasto‐plastic rheology suitable to describe the first‐order aspects of continental rock masses. Different from typical long‐term geodynamics numerical frameworks, we explicitly account for both volumetric and deviatoric response of lithospheric rocks to applied loads. Plastic correction to a viscoelastic stress state is introduced via a non‐associative Drucker‐Prager model, without resorting to the assumption of a plastic limiter. The transient behavior of crustal and lithospheric rocks is accounted for by an overstress (rate‐dependent) viscoplastic rheology, which additionally helps solving for numerical issues related to plastic strain accumulation even in the absence of energetic feedbacks. When applied to the study of the dynamics of plume‐lithosphere interactions, our implementation is able to reproduce a surface topography with complex multiharmonic wavelength patterns in agreement with observations. In the final chapter, we discuss main limitations of the current rheological description when applied to the study of transient semi‐brittle rock behavior. These aspects are tackled in a companion paper, where a thermodynamically consistent formulation extending the current numerical description is presented.


Introduction
The long-term strength of the continental lithosphere can be broadly subdivided into two end member deformation modes: pressure-dependent frictional (brittle) behavior in the upper crust following Byerlee's law; and temperature-and strain rate-dependent viscous (ductile creep) behavior in the lower crust exemplified by the Arrhenius's flow law. Assuming a (quasi)-static rheological response of lithospheric rocks to loading leads to the classical concept of the Brace and Goetze Yield Strength Envelope (YSE hereafter; Brace & Kohlstedt, 1980;Goetze & Evans, 1979). Therefore, the maximum strength of a rock is provided by the envelope given by the minimum of the two mechanisms. Following this description, with increasing temperature, that is, at greater depths, the lower crust gradually weakens upon reaching the stronger, olivine-rich lithosphere mantle. This results in a strength profile dubbed as the "jelly sandwich" YSE model. Alternatively, the crust might be able to bear all of the stress following a "crème brûlée" strength profile (Burov, 2011).
Common to both rheological descriptions is that the shape of the strength profile determines at which depth differential elastic stresses can accumulate. It therefore imposes self-consistent bounds to the amount of energy that can be stored and possibly released whether seismically or aseismically. The brittle-ductile transition at a given strain rate is given by the depth at which the two curves intersect. However, whether a rock behaves in a rather brittle or ductile fashion depends on specific environmental conditions such as the level of confining pressure, the presence and type of fluids, temperature and strain rate boundary conditions, and the total strain exerted on a rock. Experimentalists have long attempted to place some constraints on how these variables affect the rheological behavior of continental rock masses, especially under conditions spanning the transitional range between macroscopic faulting and viscous behavior (e.g. Baud et al., 2000;Brantut et al., 2013;Nicolas et al., 2016;Paterson & Wong, 2005;Pec et al., 2016;Reber et al., 2015;Reber & Pec, 2018;Violay et al., 2012;Wong & Baud, 2012;Zhu et al., 1997). The results from all these experimental 10.1029/2019JB018474 efforts led to the generalization of the concept that a certain region exists under which a rock behaves both in a ductile and in a brittle manner, or rather in a semi-brittle semi-ductile manner. Interestingly, these findings contradict the common picture of the continental lithosphere as derived on the base of the YSE concept, that is, that of strong layering as dictated by a sharp brittle-ductile transition.
To unravel the nonlinear dynamics controlling the stability, location, and extent of this semi-brittle, semi-ductile domain has become of relevance for diverse geoscientific applications. Those comprise studies aiming at describing the working of plate tectonics on Earth, subduction dynamics, and mantle-surface interactions; geological hazards finally leading to seismic ruptures; and, more recently, proper targeting of high enthalpy unconventional geothermal reservoirs.
Common to all these applications is the necessity to deal with diverse spatial and time scales over which to describe the complex nonlinear character of the underlying physics by relying on data (either from geological observations or from rock physics experiments) which are uncertain and scarce. To overcome these limiting aspects geoscientists have been relying on computer-assisted mathematical modeling techniques. Starting from the 1970s, an increasing number of sophisticated, physics-based models have been developed and applied to the study of the different Earth's system compartments, from global deformation and mantle convection (e.g. Beucher et al., 2019;Choi et al., 2013;Gerya, 2011;Heister et al., 2017;Kronbichler et al., 2012;Moresi & Lenardic, 1999;Popov & Sobolev, 2008) to basin and reservoir dynamics (e.g. Cacace & Jacquey, 2017;Cacace & Scheck-Wenderoth, 2016;Jacquey et al., 2018;Poulet et al., 2017;Räss et al., 2017;Rutqvist, 2011). However, only few attempts can be cited aiming at investigating the interactions among the different compartments and at describing the nonlinear physics occurring at their interfaces (Allison & Dunham, 2018;Barbot, 2018;Burov & Gerya, 2014;Burov & Guillou-Frottier, 2005;Cloetingh & Negendak, 2009).
As an example, the majority of models of mantle and Earth's surface interactions often neglects, or considers in a simplistic manner, the presence of the rigid upper lithosphere, which extends from the surface down to depths of (more than) 100 km. This is despite the socio-economic relevance of this geological layer, which hosts all of the world's natural and mineral resources and presents the highest hazard component.
The aim of this contribution and its companion paper is therefore to address two specific open challenges in the field of dynamic modeling applied to the Earth's system: (i) a realistic, physics-based rheological description of the lithosphere, and (ii) its flexible numerical integration. The work presented in the following will be extended in the companion study (Jacquey & Cacace, 2020) to target the multiphysics of semi-brittle, semi-ductile deformation as also relevant to strain localization and faulting processes.
In the present contribution, we describe an integrative modeling approach based on a realistic description of the lithosphere to model the full spectrum of its deformation dynamics. This includes a proper, that is, explicit incorporation of the lithosphere visco-elasto-plastic rheology including nonlinear feedback effects from the energetics of the system. Therefore, we address topics related to (i) the mathematical formulation of a visco-elasto-plastic rheology for the lithosphere; (ii) its extension to account for time-dependent brittle behavior via an overstress (viscoplastic) formulation; (iii) an implicit and efficient numerical implementation of these algorithms within a limited amount of internal iterations; and (iv) a final assessment of its numerical stability and level of performance via dedicated applications. The computational efficiency of the formulation is also quantified by a systematic analysis of the code performance against state of the art HPC environments (Jülich Supercomputing Centre, 2019) which is presented as supporting information.

Rheological Description of the Lithosphere
Depending on the pressure and temperature conditions and on the time scale of the loading process, the Earth's lithosphere can behave in a composite manner showing both a reversible (solid elastic) and irreversible (viscous and brittle) response.
In the literature, two main approaches have been put forward to describe this behavior. The first one, most commonly assumed for geodynamic problems, is based on a fluid dynamic description of the lithosphere, where the latter is approximated as an incompressible Stoke-like viscous fluid, the motion of which is controlled by a nonlinear effective viscosity. The update of the effective viscosity can include some corrections approximating brittle behavior by limiting the ductile stress resistance via an implicit plastic yielding Journal of Geophysical Research: Solid Earth 10.1029/2019JB018474 formulation, though this approach has been shown to raise some convergence issues (Spiegelman et al., 2016). The main assumption of this rheological description is that elastic stress are fully relaxed by viscous flow. This limits such an approach to modeling processes having time scales greater than the characteristic Maxwell relaxation time of the lithosphere (∼10 kyrs or higher). To avoid this limitation, modifications by integrating a deviatoric visco-elastic component have also been published (Albert et al., 2000;Gerya & Yuen, 2007;Moresi & Lenardic, 1999). While being able to account for accumulation of elastic stress, such a deviatoric formulation still lacks to incorporate a compressible volumetric component, as essential to explain dilatant and compactant rock behavior in the brittle regime (Choi et al., 2013;Yarushina & Podladchikov, 2015).
The alternative is to replace the incompressible Stoke approximation by Newton's second law of motion and to consider elasticity as the primary deformation mechanism within an explicit visco-elasto-plastic rheology (Andrews, 1978;Barbot & Fialko, 2010;Choi et al., 2013;Popov & Sobolev, 2008;Poulet et al., 2017;Yarushina & Podladchikov, 2015). This description considers the full compressibility of the material both in the deviatoric and volumetric space, where viscous and plastic corrections are applied to an initial elastic guess state. The main advantage of this approach is that it can model the full spectrum of deformation mechanisms, being elastic compressible rock behavior, nonlinear ductile creep and plastic (brittle) behavior. The latter, which is highly pressure dependent, is resolved by relying on an accurate update of the pressure considering both elastic compressibility and dilatant and compactant brittle deformation. By considering the volumetric deformation component, this rheological description allows to properly consider the nonlinearity arising from the energetics of the system (thermal feedbacks) as well as the effects of porosity and fluid pressure to the overall mechanical behavior.
When compared to a fluid approximation, a fully explicit visco-elastic-plastic rheology reveals some challenges for its numerical integration. However, it remains a valid (physically consistent) description of the lithosphere behavior under all relevant, from human to geological, time scales. All these aspects make this approach extremely flexible when applied to the nonlinear dynamics of continental lithosphere in response to both natural and anthropogenic forcing conditions.

Governing Equations
In our formulation, lithosphere deformation is approximated as a quasi-static process described by Newton's second law of motion within a continuum mechanics approximation as where i is the deviatoric stress tensor, p the pressure, the density of the material, and g the gravity acceleration vector. Equation (1) is solved for the displacement vector u as primary variable. The momentum equation is coupled to a constitutive model, which relates stress to strain and strain rate and is used to describe the rheological properties of the lithosphere.
To understand the temporal evolution of the lithosphere upon loading, we couple equation (1) to an equation describing the evolution of the energy of the lithospheric plate, therefore solving for a coupled thermomechanical problem. The energy equation, taking both heat transfer by conduction (Fourier's law) and by advection of the solid material as well as additional dissipative terms reads as where T is the temperature variable, v the solid velocity (time derivative of the displacement vector), the thermal diffusivity, H a the dissipative source term (shear heating) and adiabatic heating, respectively.
To account for the advection of material properties in a consistent way with the solid deformation history, we rely on additional compositional phase fields. Therefore, the final system is closed by solving for a conservation equation for each compositional phase as where c i is the ith compositional phase.

Constitutive Rheology
Deformation is generally described based on the velocity gradient tensor L which is obtained from the velocity v (or displacement u) vector as: The total strain rate ( . ) and the spin (w) tensors are defined as the symmetric and antisymmetric parts of the velocity gradient tensor respectively: In deriving equation (5) we assume a small strain formulation. Alternatively, the total strain rate and spin tensors can be computed based on a finite strain formulation (see Appendix A).
In this study, we consider a visco-elasto-plastic rheology. To account for accumulated strain and pressure in our rheological formulation we rely on an additive decomposition of the total strain rate tensor . . The latter is therefore expressed as a Maxwell body consisting of four main deformation mechanisms: (i) elastic, (ii) viscous, (iii) (visco)plastic, and (iv) dissipative thermal, as where the superscript e refers to the elastic component, v to the viscous component, p to the plastic component, and to the thermal component. The above described decomposition includes volumetric changes in the form of a volumetric elastic (static pressure given by the trace of the elastic strain rate tensor) and a thermal expansion component. In equation (6), each mechanism is computed according to a proper constitutive law.
Elasticity is a time-independent deformation mechanism which can be described via the Hooke's law, relating the stress tensor to the elastic strain tensor. Using a volumetric/deviatoric split as where d ∇ i dt is the Jaumann objective stress rate that satisfies the objectivity principle (does not depend on the frame of reference), which is computed using the spin tensor as The last term in the right-hand side of equation (6) represents the ability of the rock materials to dilate or compress due to a change in temperature. The expression of the thermal strain rate is therefore proportional to the rate of temperature and the volumetric thermal expansion coefficient : Ductile creep is a time-dependent deformation mechanism activated at temperature and stress levels allowing earth materials to slowly deform at constant load below their strength. Ductile creep is usually described by means of a thermally (exponential) activated (non)-Newtonian law linking the shear strain rate to the differential stress, via an effective viscosity as where e is an effective viscosity comprising diffusion and dislocation creep mechanisms:

10.1029/2019JB018474
where L and N are the diffusion and dislocation creep viscosities, respectively, anḋe L anḋe N the associated scalar strain rates and II = √ 1 2 i i is an invariant of the deviatoric stress tensor here introduced to model non-uniaxial deformation. The scalar strain rates depend on the stress and thermal states, as an Arrhenius-type of expression with laboratory-derived constants: where A L and A N are the pre-exponential constants, E L and E N the activation energies, V L and V N the activation volumes for diffusion and dislocation creep respectively. n is the power-law exponent for dislocation creep and R is the gas constant.
Lithospheric rocks have a finite strength. The build-up of viscoelastic stress due to imposed forcing can exceed the finite strength of a particular rock, which may fail in a brittle-like manner. This deformation mechanism is dominant at low temperature and low pressure. Brittle failure of rocks is often modeled using plastic deformation. Plasticity is an additional deformation mechanism which occurs when the stress state reaches a yield (strength) formulated in terms of resolved stresses.
In this study, we consider a pressure-dependent yield function following a Drucker-Prager formulation. The yield function  can be expressed as a function of the equivalent stress e = √ 3 2 i i and pressure p = − kk where = √ 3 sin ( ) and k = √ 3C cos ( ) are coefficients depending on the friction angle and cohesion C of the material. The plastic parameters and k can also be considered as linearly depending on the accumulated equivalent plastic strain to model hardening or softening of the material. We adopt a nonassociative flow rule for accumulation of plastic strain, expressed as where . is the plastic multiplier rate and  the plastic potential. The plastic potential  is formulated in terms of the dilation angle , as porous rocks generally dilate less than predicted by associative plastic models ( ≤ ): where = √ 3 sin ( ). Plastic deformation is governed by the classical Kuhn-Tucker conditions: Plasticity, whether associative or nonassociative can be used to describe time-independent brittle behavior. This description contrasts to a certain degree the outcomes of studies on kinematic friction (Barbot, 2019;Carpenter et al., 2016;Dieterich, 1978Dieterich, , 1992 and more recently experimental findings reporting a transient component of brittle deformation, even under constant stress conditions (Brantut et al., 2013;Heap et al., 2011). Additional experimental results have demonstrated dilatancy to occur at temperatures and confining pressures typical of middle to lower crustal conditions, granulite and eclogite facies (Paterson & Olgaard, 2000). All of these observations provide strong indications of an interplay between brittle-plastic and viscous behavior in rocks at far higher temperature-pressure conditions than commonly assumed. In a companion paper we will describe a more rigorous and thermodynamically consistent approach to account for these observations (Jacquey & Cacace, 2020). In the remaining of the paragraph we will only outline the fundamental concept.
In order to consider the first-order characteristics of this brittle creep behavior we augment the constitutive model by a viscous term, (overstress plasticity, Olszak & Perzyna, 1969;Perzyna, 1966). The merit of such a 10.1029/2019JB018474 formulation is two-fold. On one hand, it enables to easily account for an additional strain rate dependency during plastic deformation, as for available experimental evidences. On the other hand, even under dynamic loading conditions, the viscous term keeps the field equation hyperbolic. Therefore, strain localization in softened regions does not confine to a line of zero thickness (as usually the case for classical rate-independent plasticity), but shows a finite width in agreement with experimental data (please refer to section 4.3).
The viscoplastic strain rate is still expressed via equation (15), but the viscoplastic multiplier rate does not follow the Kuhn-Tucker conditions anymore, being calculated as where p is the plastic viscosity and ⟨x⟩ = x+|x| 2 denotes the Macaulay brackets. It is evident from equation (18) that the plastic yield can be exceeded, hence the terminology overstresses plasticity. The value of the plastic viscosity can be linked either to a physical deformation behavior or introduced only to stabilize localized deformation. The reader is referred to the companion paper for a detailed discussion on this topic (Jacquey & Cacace, 2020).
The final two equations comprising the full visco-elasto-(visco)plastic rheology of lithospheric rocks are obtained by performing a deviatoric/volumetric split of equation 6, as for the deviatoric part and for the volumetric part. Equations (19) and (20) are solved to get the stress state (p and i ) at each time step.

Numerical Implementation
In this section we outline the numerical integration of the nonlinear system of equations as described in the previous section. All the features that will be presented in the following have been integrated in an open-source dynamic package developed by the authors (LYNX, Lithosphere dYnamic Numerical toolboX).
The LYNX simulator is a Multiphysics Object-Oriented Simulation Environment (MOOSE)-based application which has been developed to model the multiphysics character of deformation with realistic lithosphere rheologies (Jacquey and Cacace, 2019).
MOOSE (Gaston et al., 2015) is an open source environment for multiphysics finite element applications developed by the Idaho National Laboratories (https://mooseframework.inl.gov/). It provides a flexible, hybrid parallel platform to solve for multiphysics and multicomponent problems in an implicit manner. It relies on state-of-the-art libraries: libMesh (Kirk et al., 2006) for capabilities related to parallel finite element method as well as geometric and I/O operations (serial and parallel I/O via ParMETIS, NetCDF, and HDF5), and a suite of scalable parallel nonlinear and linear solvers (PETSc: Balay et al., 2017;Trilinos: Heroux et al., 2005;and Hypre: Falgout et al., 2006).

Implementation of Rheology and Stress Update
To integrate the final equations for the stress-equations (19) and (20)-we follow a step-wise update procedure. This procedure is an extension of the integration presented by Dunne and Petrinic (2005) for Von Mises elastoplastic materials. The stress at the current time step is updated sequentially following an iterative procedure maintaining a precise ordering in adding the correction terms from the three rheological components, that is, first elastic, then viscous, and finally (visco)plastic.
In the following, superscript (n + 1) refers to the current time step, and (n) to the previous one.
The primary deformation mechanism accounted for in the present study is elasticity. A first, trial elastic guess for the pressure and deviatoric stress tensor is therefore computed: Journal of Geophysical Research: Solid Earth wherep e and̄e are the trial pressure and deviatoric stress tensor and̃( n) k Δt − (n) ik w k Δt the deviatoric stress tensor from the previous time step, the latter having been rotated onto the current frame.
To correct the elastic guess for viscous deformation, one needs to calculate an effective viscosity. We present in Appendix B an iterative method to update this parameter for non-Newtonian rheologies. Once having calculated the proper effective viscosity for the creep mechanisms, the equivalent viscous strain increment is obtained as 2̄e īe i the trial equivalent stress. The viscoelastic deviatoric stress tensor is then updated as follows: After correcting the trial state for the viscous deformation, a plastic correction is applied upon yielding (see equation (17)). The increment of plastic multiplier is therefore calculated by using the following identities: In order to increase the efficiency of the plastic stress update procedure we make use of the following relation: From the above equation, it follows that the plastic multiplier (for rate-independent plasticity) can therefore be calculated without the need of any iterative method as The hardening modulus H = H p + H k is related to the evolution of the plastic parameters: Δ = H Δ and Δk = H k Δ .
For viscoplastic deformation, equation (27) is modified as The final expressions of the pressure and deviatoric stress read as In LYNX we opt for iterative solution methods belonging to the class of preconditioned Krylov subspace iterative Newton solvers as provided by PETSc-SNES options. These methods are all based on successive approximations of the action of the tangent operator (also called Jacobian) on the solution vector. Through preconditioning, it is possible to improve the convergence rate of the linear system. However, the efficiency of such iterative methods still relies on a good approximation of the tangent operator from the linearized field equation. We refer to Appendix C for details on how the tangent operator for visco-elasto-plastic rheology is calculated.

Implementation of Advective Processes
Equation (3) and the advective part of equation (2) are hyperbolic scalar equations. It is well known that such equations suffer from spurious oscillations. Different numerical methods can be found in the literature to solve linear and nonlinear hyperbolic equations with high accuracy, as Riemann solvers, flux limiters, and artificial viscosity methods. The stabilization method adopted in LYNX is based on the Entropy Viscosity formulation as presented by Guermond et al. (2011). Our choices stem from the fact that this method is relatively easy to implement in the working framework, is independent of the mesh type and order of polynomial discretization used, and the dissipative viscosity term introduced is consistent with the principle of minimum entropy.
The basic idea of the method is to stabilize the equation by adding a nonlinear dissipative term in the form of a viscous flux (g , where u h is the numerical approximation of the exact solution of the advected quantity at the current time) that satisfies the minimum entropy principle ( ≤ 0). The artificial viscosity ( ⋆ ) is therefore obtained on the base on the local size of the entropy production, as proportional to the absolute value of the entropy residual limited by an additional viscosity proportional to the local maximum eigenvalue of the stiffness matrix.
For both the temperature and compositional phases, the temporal integration is based on a second-order, A-stable Backward Differentiation Formula. To insure stability of the Entropy Viscosity method (Guermond et al., 2011), we impose an upper threshold for the value of the time step following a classical Courant-Friedrich criterion. An example of the advection scheme is presented in Figure 1. An animation of this benchmark test is also provided as supporting information.

Results
In this section we present a set of synthetic and realistic applications intended to test the ability of the simulator to deal both with single-and multiple-deformation mechanisms. In doing so, we first present two applications for which we could derive analytical solutions. These will serve to benchmark the simulator. In a final stage, we describe some realistic applications, including localization under shear and dynamic coupling between plume dynamics and surface deformation within a rheological consistent lithosphere.

Stress Build-up in a Maxwell Visco-elastic Medium
To validate the implementation of the visco-elastic stress update, we consider the following benchmark, as originally introduced by Gerya and Yuen (2007). We consider a Maxwell visco-elastic medium (shear modulus G = 1.0 × 10 10 Pa and viscosity = 1.0 × 10 22 Pa s) under pure shear at constant strain rate ( . = 1.0 × 10 −14 s −1 ). The analytical solution to the problem is given by  The agreement between the numerical and analytical results validates the implementation of the stress update procedure for visco-elastic rheology and demonstrates the accuracy of the numerical results.

Uniaxial Elasto-(Visco)Plastic Loading
This benchmark is intended to verify the implementation of the plastic update. We consider an elasto-(visco)plastic Von Mises medium (G = 2.0 × 10 9 , K = 2.0 × 10 9 , and C = 1.0 × 10 7 Pa) subjected to uniaxial loading at constant velocity (v = 1.0 × 10 −5 m/s). Following the setup, the deviatoric stress builds up with time until it reaches the finite strength of the material, which then yields plastically. We consider three cases here, one relying on perfect plasticity (no hardening), another with linear hardening (hardening modulus of H = 5.0 × 10 8 Pa), and a last one with viscoplastic deformation with hardening (plastic viscosity p = 4.0 × 10 11 Pa s for the same hardening modulus as before). For the three cases, we derive the analytical solutions for the equivalent stress ( e = √ 3 2 i i ) as • Perfect plasticity: • Plasticity with constant hardening: • Viscoplasticity with constant hardening (t p = √ 3C 3Gv corresponds to the onset of plastic deformation): Figure 3 shows the model setup as well as the evolution of the equivalent stress for the three cases. At the beginning, the equivalent stress increases linearly as elastic deformation is accumulated within the domain. When the equivalent stress reaches the cohesion of the material, plastic deformation starts to accumulate. For perfect plasticity, the stress cannot exceed the constant strength of the material and therefore stays constant. Including hardening, the strength linearly increases with the accumulated plastic strain, as the equivalent stress. When accounting for viscoplastic deformation in addition, the equivalent stress exhibits some nonlinear hardening before increasing linearly with the accumulated viscoplastic strain. For the three cases, the numerical solution matches exactly the analytical solution, demonstrating the robustness of the numerical implementation of the (visco)plastic update.

Shear Bands in a Visco-Elasto-(Visco)Plastic Medium
Shear bands are localized regions of high strains forming by increased inelastic flow. Their formation is likely related to the presence of some material heterogeneities (weak domains). Due to increased strain rates, shear bands are important precursors to macroscopic failure events. Therefore, much attention has been given to the study of onset conditions of shear bands and their evolution. While linear stability analysis can be effectively used to describe the stability of such deformation regions, it fails in describing the nonlinear physics behind their evolution. An alternative is to rely on numerical modeling techniques.
Generally, ductile shear zones are modeled by considering thermal softening of the material by coupling the energetic of the bands to their viscous deformation (last term in equation (2)). Additional effects from grain size reduction might also contribute to the softening. The energy dissipated by the irreversible deformation causes the temperature to increase, which results in the build-up of strain rates, power law decrease in the effective creep viscosity, and further softening and localization. The formation of a viscous band follows two stages, that is, an initial narrowing phase due to inelastic strain increase, and a subsequent widening of the band by thermal diffusion and viscous heating (see Figure 4; Duretz et al., 2014). Therefore, considering shear bands as thermo-mechanical viscous instabilities requires strain softening by thermal dissipative The geometry and boundary conditions used for this example are described in Figure 5a. processes for localization to be maintained. However, power law viscous rheology or grain size-induced softening under linear viscous flow might not always lead to effective localization over time. In addition, as demonstrated by Regenauer-Lieb and Yuen (2003) the onset of strain localization and shear bands formation can be consistent with strain hardening. An additional limitation of a purely viscous approach to the formation of shear bands relates to the viscosity contrast between the highly deformed band and the surrounding matrix. Such a gradient will induce overpressure within the band, which will favor fluid to escape into the matrix contrasting the presence of channelized fluid as commonly observed within mylonitic shear zones (Bauer et al., 2000;Spruzeniece & Piazolo, 2015;White et al., 1980).
To compensate for this aspect, viscous flow rheology can be augmented with a pressure dependent plastic component in the form of classical Drucker-Prager or Mohr-Coulomb rate-independent formulation. In this example, we consider a setup introduced by Duretz et al. (2018) consisting of a 4 by 2 km domain subject to pure shear at constant strain rate. A circular weak inclusion is considered at the center of the domain to perturb the stress field and initiate strain localization. The model geometry and boundary conditions are illustrated in Figure 5a). Whether a viscoplastic approach to ductile shear zones formation might solve for the pressure gradient inside the band, it also comes with some problematics. The major challenge while relying on a classical continuum mechanics approach, stems from the lack of any internal length scale of the process, that is, the instability localizes into a domain of zero thickness. Indeed, rate-independent plastic formulations require the bands to narrow indefinitely in contrast to the majority of observations (see Figure 5b). Overstress viscoplasticity by definition introduce a characteristic internal length scale and therefore can prevent the strain from localizing into infinitely narrow shear bands. For such models, the width of the shear zones will not only depend on the size of the imperfection initializing the localization but also on the internal characteristic length, which has an inverse dependence on the viscosity parameter ( p in Figure 6. Model geometry for the plume-lithosphere example. equation (18)). This is illustrated in Figure 5c and Figure 5d. When the size of the imperfection is smaller or of the same order of the internal length scale (Figure 5c), the width of the resulting band is dominated by the former parameter, leading to a final results that resembles that obtained with rate-independent plasticity. By increasing the plastic viscosity (see properties in Table 1), therefore decreasing the internal material length scale, the amount of accumulated plastic strain is decreased and the amount of overstress is increased. This leads to a diffusion of the localized deformation and therefore to a widening of the shear bands.

Influence of a Rising Plume on the Lithosphere Deformation
The last example deals with modeling the interaction between deep mantle processes and lithosphere deformation. At this purpose, we consider a setting that has been extensively investigated within the geodynamic modeling community, but for which no general consensus has been achieved (Burov & Gerya, 2014;Burov & Guillou-Frottier, 2005;d'Acremont et al., 2003;François et al., 2017;Guillou-Frottier et al., 2007;Koptev et al., 2015Koptev et al., , 2016. The problem is that of a rising plume and its thermo-mechanical coupling to the deformation of a visco-elasto-plastic lithosphere. The plume is modeled as a density anomaly having both a thermal and a chemical signature. By combining a complex rheological description, thermo-mechanical coupling and given the large scale of the model, this study is a good proxy to finally evaluate the performance of the modeling approach presented in this contribution. The geometry of the model is presented in Figure 6. We consider a two-dimensional domain of size 2,000 km by 650 km which comprises a 40 km crust (20 km upper crust and 20 km lower crust), a 100-km-thick lithospheric mantle and an underlying asthenospheric mantle. Consistent to Burov and Guillou-Frottier (2005), the plume is represented by a spherical thermal and density anomaly of 300 km diameter centered at the bottom of the model.
The deformation of the different geological layers is described by a visco-elasto-plastic rheology with thermally activated non-Newtonian viscosities. The parameters adopted for this example are summarized in Table 2. Creep parameters assume quartz, diabase, and olivine as main rheologies for the upper crust, lower crust, and lithospheric mantle, respectively. The plume rheology is similar to that of the mantle, but differs in terms of viscosity because of the higher temperature condition. We also apply a free surface boundary condition at the top of the model, in order to be able to track the evolution of the surface topography as a Creep -activation energy E N 5.20 × 10 5 J/mol result of the plume-lithosphere interactions. In order to maintain the current formulation consistent with previously published studies (Burov & Guillou-Frottier, 2005;d'Acremont et al., 2003), in this example we neglect any heat production term in the right-end side of equation (2).
A lithostatic pressure distribution is assumed as initial conditions, based on the density distribution. The initial temperature is derived from a solution to the half-space cooling model (d'Acremont et al., 2003), considering a 50 Myr old lithosphere: where is the depth within the lithosphere, T b = 1,330 • C is the temperature at the bottom of the lithosphere, is the thermal diffusivity, and t is the age of the lithosphere. In the asthenospheric mantle, we adopt a near-adiabatic temperature gradient up to 2,000 • C at the base of the model.  Figure 7 illustrates the initial temperature profile as well as the resulting strength profile, assuming a background strain rate of 1.0 × 10 −15 s −1 . The initial temperature within the plume is increased by ΔT = 250 • C being also characterized by a density anomaly of Δ = Δ + Δ c = 25.5 + 25 = 50.5 kg/m 3 due to thermal and chemical anomalies, respectively (Burov & Guillou-Frottier, 2005).
In Figure 8 we illustrate the evolution during the first 10 Myr in terms of the distribution of the compositional phases, temperature, and effective viscosity, the latter taken as a proxy for the strength of the different layers.
The density anomaly in the plume triggers buoyancy forces leading to its rising. Deformation within the lithosphere results from the combination of two driving forces: (i) build-up of elastic stress in response of the rising of the plume, and (ii) their relaxation by viscous flow activated by the thermal anomaly associated with the plume itself.
As the plume head reaches the lithosphere-asthenosphere boundary, the mantle lithosphere thins in response to compressional stress. The amount of thinning is up to 50% of the initial thickness in the center of the model domain beneath the plume head. However, compression in the mantle extends laterally over a larger area than the plume head itself, around 600 km. During these early stages in the evolution of the plume-lithosphere system, the lithosphere still maintains a finite strength as indicated by the resulting viscosity map. Within a relative strong lithosphere, stress cannot be relaxed by viscous flow, but rather accumulate (visco-)elastically. This finally leads to flexural bending of the lithosphere and to a biparted mechanical response of the plate, with compression around a limited domain centered on the plume head and extension at its margins. The arrival of the plume at the lithosphere-asthenosphere boundary, at around 1 Myr, marks the onset of lithosphere weakening. Strain starts to accumulate and local temperature conditions to increase. This thermo-mechanical destabilization leads to a decrease in the effective viscosity, though at this stage, limited to a domain centered around the plume head. Despite the observed weakening, the lithosphere still maintains some strength, as indicated by ongoing build-up of stress at the plume's edges. As the plume starts to flatten across the lithosphere-asthenosphere boundary, the mechanical response of the overlying lithosphere changes. Flattening of the plume finally results in a mechanical decoupling at middle crustal levels. The rheological configuration of the lithosphere resembles that typical of a "jelly-sandwich" lid, with a relative strong upper crust and mantle, bearing stress visco-elastically, decoupled by the presence of a weak middle crust.
The continental lithosphere acts as a mechanical buffer to the plume dynamics. Therefore, variations in the mechanical configuration of the plate reflects changes in the surface expression. These are exemplified in Figure 9, where the evolution of the surface topography resulting from these interactions are plotted.
The first surface expression of the plume-lithosphere interaction is a major mantle domal uplift reaching a maximum topographic elevation of 1 km. This relative uniform uplift has a primary wavelength of 800 km. Due to the onset of weakening in the lower-to-middle crust, there is an attenuation of the principal wavelength, and onset of subsidence localized to a domain of finite extent centered around the plume. The region characterized by downward surface deformation is flanked by adjacent domains where the surface topography reaches a maximum height. This trend of upward and downward surface movements is a feature that characterizes the whole of the evolution of the system, with the amount of subsidence gradually increasing with time. Horizontal flattening of the plume and corresponding weakening of the overlying lithosphere, finally lead to a mechanical decoupling at middle crustal levels, and result in a further attenuation of the principal wavelength. The heterogeneous response of a mechanically decoupled lithosphere to the stress induced by the plume impinging at its base leads to the superposition of secondary wavelengths (30 to 50 km and 80 to 100 km) onto the first characteristic. The final topography, therefore, does not show a single, large-scale, and monoharmonic uplifted area, but has a more heterogeneous signature consisting in distributed basins of different characteristic wavelengths, which more accurately resembles available observations (

Discussion
In this contribution, we have described an explicit implementation of a visco-elasto-(visco)plastic rheology for modeling deformation of the continental lithosphere. In the current formulation, brittle deformation is accounted for via a (visco)plastic stress correction including linear hardening/softening of the yield function.
A proper understanding of the physics behind rock deformation under crustal temperature and pressure conditions is an essential component for any geodynamic numerical toolbox. This, in turn, requires to overcome assumptions that are usually done while describing the response of rocks under dynamic loading in order to reconcile the amount of experimental and field evidence that is nowadays available. What follows is a brief discussion on what we consider as major open questions in this respect. Emphasis will be on limitations of current geodynamic numerical tools when applied to the description of brittle rock behavior.
In the companion paper, we will present our efforts to extend the current formulation to overcome these assumptions and limitations (Jacquey & Cacace, 2020).
Brittle deformation is usually described as a time-independent process, approximated by a particular flavor of rate-independent plasticity, J2-type or Mohr-Coulomb. However, an increasing number of laboratory experiments have been starting to report an intrinsic transient response of (drained and undrained) brittle rocks (Blanpied et al., 1995;Kilgore et al., 1993;Toro et al., 2004), even under constant stress conditions (Brantut et al., 2013;Heap et al., 2011). These transient effects are thought to influence subcritical crack growth, the precursory dynamics to earthquake rupture and volcanic eruptions, the effective recovery of geothermal reservoirs, and the safety of long-term storage of hazardous materials in the underground. Therefore, a proper quantification of the causative physics of time-dependent (semi) brittle behavior has not only a scientific merit but also a socio-economic relevance.
The complex behavior of brittle rock as reported by the majority of triaxial deformation experiments is difficult to reconcile by relying on empirical laws, where hardening in the rock cohesion and friction properties are depicted as linear functions of the resolved accumulated plastic strain. Though the physical grounds of the resulting dynamics is still a topic of ongoing research, some degree of consensus begins to behold in linking this behavior to first-order processes related to the evolution of microstructural discontinuities and defects, the type of fluid permeating the pore space and fluid-mediated reactions with the rock skeleton.
Laboratory experiments on rock samples under isotropic compression have reported inelastic strain accumulation accompanied by a clustering of acoustic emissions and an acceleration in their frequency spectra upon reaching a critical pressure. The value of this critical pressure for inelastic strain accumulation has also been shown to structurally depend on the initial nominal porosity of the rock investigated and its evolution during the loading experiment. Based on these observations, previous studies have concluded that porous rocks exhibit a finite volumetric strength and that phenomena as pore collapse or grain crushing are intrinsically linked to the evolution of the rock's pore network. However, normally adopted plastic yield functions assume the strength of a rock to be a linear function of the lithostatic pressure only.
In the companion paper, we will extend the current rheological description to include processes that can reconcile the aforementioned observations. Based on thermodynamics considerations, we will show how to derive a physical formulation to incorporate processes that are usually neglected in geodynamic modeling, including (i) rock weakening and localization of deformation via the introduction of a consistent damage rheology, (ii) the additional impact of pore pressure on the bulk deformation for a two-phase fluid-rock matrix system, (iii) the role of porosity and its evolution in the constitutive mechanical behavior of porous rocks, and (iv) when coupled to the damage history of microdefects, its precursory role to macroscopic failure.

Conclusion
In this contribution, we have introduced a novel simulator for modeling thermo-mechanically coupled processes driving the deformation dynamics of the lithosphere. The LYNX simulator relies on the open-source massively parallel finite element-based numerical framework MOOSE to solve for the nonlinearities of the coupling among the different processes and inherent to the visco-elasto-(visco)plastic rheology adopted. The rheology is implemented relying an explicit constitutive update with a limited amount (only for non-Newtonian viscosity) of iterations. We illustrated the capability, robustness, and parallel performance of the code using a suite of numerical examples of increasing complexity, ranging from analytical benchmarks to crustal scale example.
Improving our understanding of the driving forces in long-term tectonic models requires a better description of the individual deformation mechanisms, their nonlinear coupling and of the resulting rheology. Elastic stress accumulation and compressibility are essential processes of the lithosphere motion and the compactant/dilatant behavior during plastic deformation can significantly alter the bulk pressure response. In this regard, the current modeling framework provides a powerful tool to analyze the dynamic effects of brittle and ductile mechanisms on the overall behavior of the lithosphere.
In the following we describe how to modify the formulation presented in the body of the text to account for a finite strain formulation. Finite strain is based on the definition of the deformation gradient F, defined as: where I is the identity tensor and u the displacement vector.
The velocity gradient (L) is then calculated using The strain rate and spin tensors are then calculate as described by equation (5). Using a finite strain formulation allows to account for second-order terms in the strain update which are neglected in a small strain approximation.

Appendix B: Iterative Update of the Effective Viscosity
When using a non-Newtonian viscosity (dislocation creep), the viscous strain increment depends on the current stress conditions. It is therefore appropriate, if not required to iterate the update in order to find a stable solution. In this context, assuming a dislocation creep rheology, we aim at finding the zero of the following residual: At the beginning of the update procedure, Δe v = 0 is determined based on the elastic trial guess. The algorithm then searches for a root of the above equation. Several algorithms deal with root-finding of single-variable function. In LYNX, we relied on two algorithms as described in (Press et al., 2007), that is, (i) Brent's method (Brent, 1971), and (ii) a modified Newton-Raphson procedure. For the latter, a derivative of the residual is required, which is expressed as

Appendix C: Consistent Tangent Operator
An important component of the simulator to insure numerical convergence is an appropriate definition of the consistent tangent operator for the visco-elasto-(visco)plastic rheology described in this contribution. Using the stress expression defined in equations (29) and (30), the consistent tangent operator is defined as where D ve is the viscoelastic consistent tangent operator, I i kl = 1 2 ( ik l + il k ) is the symmetric identity tensor and E p is the plastic correction to the consistent tangent operator.
The viscoelastic consistent tangent operator is expressed as where I vol i kl = i kl and I dev i kl = I i kl − 1 3 I vol i kl .

Journal of Geophysical Research: Solid Earth
where vp is the viscoplastic correction as defined in equation (28).