Multiphysics Modeling of a Brittle‐Ductile Lithosphere: 2. Semi‐brittle, Semi‐ductile Deformation and Damage Rheology

The brittle‐ductile transition is a domain of finite extent characterized by high differential stress where both brittle and ductile deformation are likely to occur. Understanding its depth location, extent, and stability through time is of relevance for diverse applications including subduction dynamics, mantle‐surface interactions, and, more recently, proper targeting of high‐enthalpy unconventional geothermal resources, where local thermal conditions may activate ductile creep at shallower depths than expected. In this contribution, we describe a thermodynamically consistent physical framework and its numerical implementation, therefore extending the formulation of the companion paper Jacquey and Cacace (2020, https://doi.org/10.1029/2019JB018474) to model thermo‐hydro‐mechanical coupled processes responsible for the occurrence of transitional semi‐brittle, semi‐ductile behavior in porous rocks. We make use of a damage rheology to account for the macroscopic effects of microstructural processes leading to brittle‐like material weakening and of a rate‐dependent plastic model to account for ductile material behavior. Our formulation additionally considers the role of porosity and its evolution during loading in controlling the volumetric mechanical response of a stressed rock. By means of dedicated applications, we discuss how our damage poro‐visco‐elasto‐viscoplastic rheology can effectively reconcile the style of localized deformation under different confining pressure conditions as well as the bulk macroscopic material response as recorded by laboratory experiments under full triaxial conditions.


Introduction
At a macroscopic scale, irreversible rock deformation mechanisms are subdivided into either brittle or ductile, depending on the capability of a rock specimen to accommodate substantial strain with or without developing macroscopic fractures. The failure mode of both porous and compactant rocks undergoes a transition from brittle fracturing to ductile flow upon increasing confining pressure and/or ambient temperature. Indeed, the brittle-ductile transition is usually considered to represent a region at depths where the deformation of lithospheric rocks switches from a brittle-pressure-dependent frictional-behavior to a ductile-temperature-and strain-dependent-behavior. Experimental advance accumulated in the last decades have demonstrated that such transition is not as sharp and stable as dictated by the classical Yield Strength Envelope model of Goetze and Evans (1979) and Brace and Kohlstedt (1980), being the result of a rather complex interaction among different microscopic processes sensitive to the specific external conditions (i.e., confining pressure, temperature, strain and strain rates, and presence of fluid phases) (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;Wong & Baud, 2012;Violay et al., 2012;Zhu et al., 1997). It follows that the brittle-ductile transition in porous and compactant rocks comprises a finite region where the rock gradually transitions between the two end-member deformation modes showing a semi-brittle, semi-ductile response to applied loads.
Unraveling 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. Developing a physical framework relevant to these applications requires accounting for diverse spatial and time scales over which to describe the coupling and non-linearities between the different physical processes while relying on observations, whether from laboratory experiments or geological observations. In order to describe the physics behind these observations, different modeling approaches have been developed over the last decades, including Cosserat continuum Rattez, Stefanou, Sulem, Veveakis, et al., 2018), different flavors of capped plastic models (Grueschow & Rudnicki, 2005;Poulet & Veveakis, 2016;Rutter & Glover, 2012), phase field approaches (Bryant & Sun, 2018;Sun et al., 2011), a poro-visco-elasto-plastic model including a bulk viscosity by Yarushina and Podladchikov (2015), and damage rheologies (Hamiel, Liu, et al., 2004;Karrech et al., 2011Karrech et al., , 2011Lemaitre, 1985;Lyakhovsky et al., 1993). Common to all these efforts is to arrive at a consistent integration of micromechanical aspects into a macroscale averaged constitutive model of the rock sample. Mechanical models that rely on a damage rheology aim at accounting for the evolution of grain sizes or contact areas on the macroscopic scale by degrading the elastic moduli and accelerating inelastic accumulation.
In the present contribution, we extend the physical framework presented in Jacquey and Cacace (2020) to target the multiphysics of semi-brittle, semi-ductile deformation as also relevant to localization of deformation and faulting processes. Therefore, we describe how to couple the microstructural evolution, expressed in terms of a damage intensity variable, onto the macroscopic deformation and frictional behavior of porous rocks by also accounting for the presence of fluids via the introduction of a dynamic porosity. The modeling results will be finally used to demonstrate how fluid-assisted damage weakening intrinsically leads to macroscopic localization and how these non-linear feedbacks might control dilatant brittle-like deformation.

Mathematical Formulation of the Problem
The mathematical formulation presented in this contribution is based on and further expands the explicit visco-elasto-viscoplastic rheology described in the companion paper Jacquey and Cacace (2020). The evolution of damage in the rock upon varying loading conditions is based on the work of Hamiel, Liu, et al. (2004), , and Lyakhovsky and Hamiel (2007) which enables to include the effects of fluid flow in porous rocks, thus resulting in a tightly coupled thermo-hydro-mechanical problem.

Governing Equations
In this section, we describe how the governing equations for thermo-mechanical coupling as introduced in Jacquey and Cacace (2020) can be reworked to account for hydro-mechanical coupling and damage rheology. The final coupled system of equations is shown to describe the deformation of rock materials, their degradation due to progressive damage of the rock framework as also coupled to the energetics of the two-phase pore fluid and rock matrix system.
The deformation of porous rock materials 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 effective pressure, p the fluid pressure, B the Biot poroelastic coefficient, the density of the material, and g the gravity acceleration vector.
In equation (1), we make use of the effective stress law introduced by Carroll and Katsube (1983), consistent with the Biot's poroelastic theory (Biot (1941(Biot ( , 1973), which relates the effective stress tensor σ ′ , the Cauchy stress tensor σ, and the fluid pressure p as σ ′ i = σ i + B p i . The effective pressure p ′ is therefore calculated The presence of pores within a rock allows fluids to circulate and to interact with the solid matrix. In this study, we consider the pores to be fully saturated by a monophasic fluid (liquid water in the most general case). Fluid flow within the pores of geomaterials is generally described by means of Darcy's law, which is derived from the Stokes equation by means of homogenization over a Reference Elementary Volume. To derive the governing equation for the fluid pressure, we enforce the mass balance for both the fluid and the solid phases as: where is the porosity of the material, v and v s the fluid and solid velocities, and and s the fluid and solid densities. From the fluid mass balance (equation (2)), one can derive the governing equation for the fluid pressure: where M B is the Biot modulus or reciprocal of the storage coefficient, e kk and in kk the elastic and inelastic volumetric strain, respectively, and q D the Darcy's velocity or specific discharge. The latter is derived from the conservation of the linear fluid momentum leading to the following definition: where k is the permeability tensor and the fluid viscosity.
The solid mass balance (equation (3)) can be used to cast the governing equation for the evolution of the porosity as 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 i th compositional phase.

Thermodynamics for Damage Rheology
The damage rheology adopted in this study is based on thermodynamics considerations as introduced in Lyakhovsky et al. (1993) and . In what follows we describe only the main points behind the formulation, while the interested reader is referred to the literature for a more comprehensive explanation of its derivation.
We consider the following thermodynamic state variables describing the state of a damageable, porous, and fully filled rock: (i) the Cauchy elastic strain tensor i , (ii) the damage intensity variable , (iii) the variation of fluid volume content , and (iv) the inelastic volumetric strain in v = in kk . The free energy of a porous and fully saturated damaged rock can be expressed as a sum of the elastic energy under drained conditions Ψ dr and of the poroelastic coupling term of the saturated medium Ψ (Hamiel, Liu, et al., 2004;:

Journal of Geophysical Research: Solid Earth
where e v = e kk and e d = √ 2 3 e e i e e i are the volumetric and deviatoric elastic strain invariants, || e i || = √ e i ∶ e i is the norm of the elastic strain tensor, K ( ), G ( ) are the damage-dependent bulk and shear moduli, and Γ ( ) is a third elastic modulus introduced to describe the evolution of the stiffness of damaged materials. The deviatoric elastic strain tensor e e i is expressed as e e i = e i − 1 3 e v i . In a similar way as discussed by Lyakhovsky et al. (2015), we adopt the following description of the elastic moduli as a function of damage: where 0 is the critical elastic strain ratio or modified internal friction, considered here as a material constant.
The Cauchy stress tensor σ i and the fluid pressure p are defined as thermodynamic forces associated with the elastic strain tensor e i and the volume fluid content , respectively: where in equation (11) we recovered the relation between the Cauchy stress tensor, the effective stress tensor, and the fluid pressure. The equivalence of equations (4) and (12) for evolving fluid pressure is demonstrated in Appendix B. The effective stress tensor σ ′ can therefore be obtained based on equations (11) and (12) as whereσ ′ is the undamaged effective stress and σ the damaged stress correction, expressed as where = e v || e || is the elastic strain ratio and is bounded from − √ 3 for isotropic compression to + √ 3 for isotropic tension.

Flow Laws and Damage Rheology
In a similar way as discussed in the companion paper, also in this study we consider the porous rock as a Maxwell body consisting of (i) an elastic component, (ii) a viscous correction, (iii) a viscoplastic (in an overstress sense) term, and (iv) a final dissipative component related to the energetics of the rock. The total strain rate tensor can therefore be split as where the superscripts e refers to the elastic part, v to the viscous part, vp the viscoplastic (overstress) part, and the thermal part of the strain rate tensor.
The viscoplastic correction to the viscoelastic guess state is considered in the form of an overstress plastic model in the undamaged state (σ). The yield function is formulated in terms of the projections of the stress invariants onto the yield surface as 10.1029/2019JB018475 wherep ′0 and̄0 e are the projections of the undamaged effective pressure and equivalent stress, respectively, on the yield function. Note that from equation (17), the yield function is formulated in the undamaged stress space, which is equivalent to the elastic strain formulation as adopted by Lyakhovsky et al. (2015).
The viscoplastic strain rate is expressed as described in Jacquey and Cacace (2020): where p is the plastic viscosity and ⟨x⟩ = x+|x| 2 denotes the Macaulay brackets. Using the general expression of the yield function introduced in equation (17), one can express the volumetric and deviatoric viscoplastic strain rates as Damage accumulation occurs when viscoplastic strain accumulates and is therefore expressed as an overstress quantity by where is the damage viscosity.
In this study, we distinguish between the plastic and damage viscosities, each being considered as an independent material property. Our choice stems from the fact that the functional dependency of the two viscosities upon variations in the internal primary variables (pressure and temperature) is qualitatively different. Therefore, in our formulation, the damage intensity variable cannot be directly associated to the accumulated viscoplastic strain. The main difference between the formulation presented above and the one from Hamiel, Liu, et al. (2004) and Lyakhovsky et al. (2015) is that we adopt a rate-dependent plastic model (overstress plasticity) instead of their viscoelastic formulation. The reasons behind our choice are related to the fact that by our formulation we can (i) easily capture the full dependency on loading rates as compared to classical rate-independent plastic models by (ii) keeping an explicit expression of the strength of the material as compared to a viscoelastic model.

Onset of Inelastic Deformation and Damage Accumulation 2.4.1. Low-Porosity Brittle Materials
The inelastic behavior of low-porosity materials has been historically described by relying on classical Mohr-Coulomb failure analysis. The latter is formulated in terms of a functional dependence of the yield upon the friction coefficient to include a pressure-sensitive correction to the strength of rock materials. For low-porosity materials (nominal porosity less than 3% to 5%), friction coefficients are defined to follow Byerlee's law (static friction). More advanced formulations include additional dynamic corrections based on empirical weakening/strengthening laws (e.g., rate-and-state formulations). In the present contribution, we rely on a static friction definition while hardening and weakening effects are accounted via the overstress plastic model and the damage rheology, respectively. We describe the plastic yield function for low-porosity materials in terms of the elastic strain ratio as introduced by Lyakhovsky et al. (1997), as The plastic yield can therefore be expressed in terms of the undamaged stressσ in the form of a cohesionless Drucker-Prager plastic yield formulation as The critical elastic strain ratio is usually referred to as internal friction coefficient of the material and can be evaluated based on an apparent friction angle (Lyakhovsky et al., 1997). It follows that with this formulation we can directly compare the friction coefficient from classical approaches to the critical elastic strain ratio 0 adopted here.

High-Porosity Semibrittle, Semiductile Materials
The mechanical behavior of porous rocks is controlled to a first degree by their initial porosity and its evolution due to the changing stress conditions. Porous sandstones exhibit evidences of inelastic accumulation-non-linear stress-strain relation accompanied by emission of acoustic emissions-also during isotropic compression tests (Wong & Baud, 2012). These observations provide more than indications on the existence of a (porosity-dependent) critical pressure capping the elastic domain.
We model such behavior by making use of a capped yield function which depends on the actual porosity in the form (Hamiel, Liu, et al., 2004;Lyakhovsky & Hamiel, 2007) where D ( ) is a porosity-dependent coefficient. The critical strain ratio can also be expressed in terms of the undamaged stress invariants as where the critical effective pressurep is the (undamaged) pressure at which inelastic deformation occurs under isotropic compression. The yield function can therefore be formulated in a similar way to equation (22) by substituting 0 with cr .
A major problem with the yield function described above is that it is not defined in certain regions of the undamaged stress space. This is due to the fact that the critical strain ratio introduced in equation (24) decreases when confining pressure increases. To overcome this issue, we reformulate the yield function in a more general sense. Following the new formulation, the yield function is considered to depend on the values of the projected stress onto the same yield, the latter being calculated based on the square of the yield function itself. A more comprehensive description on how this procedure has been implemented is detailed in section 3.2. Figure 1 illustrates the yield function for (a) low-porosity-frictional-materials and (b) high porous rocks in the undamaged stress space. As specified by Lyakhovsky et al. (2015), adopting a undamaged stress formulation or elastic strain formulation for the yield function of materials allows to account for material weakening without the necessity of updating the yield function as it is usually done if relying on a stress space formulation.

Numerical Implementation
In this section we outline the numerical implementation adopted to solve for the non-linear system of equations as described in the previous section. All of the features that will be presented in the following have been already integrated in the open-source dynamic package LYNX (see also the companion paper at this purpose). LYNX (Lithosphere dYnamic Numerical toolboX) is a MOOSE-based application developed to model the multiphysics character of deformation based on a realistic description of the rheological behavior of the lithosphere Jacquey and Cacace (2019). MOOSE (Multiphysics Object-Oriented Simulation Environment) (Gaston et al. (2015)) is an open-source environment for multiphysics Finite Element applications developed at 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).
In the next paragraphs, we detail the update procedure for stress and elastic strain when considering accumulation of damage and the integration algorithm for the convex implicit yield formulation for high-porosity materials.

Update Procedure for Damaged Materials
The damage rheology adopted in this study is formulated in terms of elastic strain. An important aspect of our approach is that the elastic moduli structurally depend on the damage intensity variable. Therefore, the current formulation can be easily adopted to account for stress as primary rheological parameter.
The undamaged elastic strain at the current time step is updated sequentially following an iterative procedure maintaining a fixed order in adding the specific correction terms from the three rheological components, that is, first elastic, then viscous and finally (visco)plastic. Only at the final update step, the stress tensor is modified based on the elastic strain tensor and the damage-dependent stiffness. In the following, superscript (n + 1) refers to the current time step, and (n) to the previous one.
Following the description in the companion paper Jacquey and Cacace (2020), the primary deformation mechanism accounted for is elasticity. The trial elastic strain is therefore computed as wherẽe (n) i is the elastic strain from the previous time step rotated to the current frame and Δ (n+1) is the total strain tensor at the current time step. For the viscoelastic update of the material, the reader is referred to Jacquey and Cacace (2020).
To compute the viscoplastic update, we calculate the undamaged stress invariants based on the trial (visco)elastic strain asp The viscoplastic strain increments are derived from the two invariants by integration as are the viscoplastic correction for the volumetric and deviatoric parts, respectively. The resulting elastic strain increment is then updated using the following relationship: Finally, the stress tensor is updated using a Taylor expansion of equation (13) with respect to the elastic strain tensor and damage variable as Figure 2. Update of the reference point or pressure to evaluate the projections of the stress invariants on the yield surface. The reference pressure corresponds to the dilation-compaction limit based on the given trial deviatoric stress. In the figure, the updated reference pressure is shown for two end-member cases, that is, in the dilation and compaction domains under the same trial deviatoric stress.
whereσ ′(n) i is the effective stress tensor from the previous time step rotated onto the current frame and D i kl defined as where C i kl is the elastic stiffness evaluated for the undamaged moduli and 1 i kl is the symmetric identity tensor.

Evaluation of a Convex Implicit Yield Function
In section 2, we have introduced the flow laws to describe both inelastic strain and damage accumulation. Both laws depend on the evaluation of the yield function. The latter has been formulated in terms of the projected undamaged stress invariants onto the yield surface. If based on simple yield surfaces (Mohr-Coulomb or Drucker-Prager), the stress invariants projection can be computed explicitly following similar working steps as for updating the stress invariants under a perfectly plastic model. Therefore, one can follow the procedure as described in Jacquey and Cacace (2020). However, when dealing with a convex yield surface, such as the one described in section 2.4, it is no longer straightforward to calculate the projected quantities. Therefore, one must rely on specific algorithms.
In this study, we rely on a procedure modified after Stupkiewicz et al. (2014), which makes use of a reference point located inside the convex yield surface. The reference point defines an elastic state and can be used together with the trial state (in the viscoplastic region) to find the projections of the stress invariants. We define this point in relations to the trial state, which makes our modified procedure equivalent to a non-associative version of the original approach described in Stupkiewicz et al. (2014). The reference point is here defined by assuming zero deviatoric stress and a pressure value corresponding to the transition from dilation to compaction under the given trial deviatoric stress. In case the computed reference pressure if found to exceed the boundaries of the capped yield, we make use of the critical pressure instead. Figure 2 schematically summarizes the geometrical update procedure.
Once the reference point has been found, one can evaluate the invariants projections onto the yield surface. From Figure 2, it is clear that the projection point in the stress invariants space is aligned with the trial state and the reference point. By introducing as the distance from the reference point, we have the following notations following Stupkiewicz et al. (2014):p wherep ′ r is the reference effective pressure and e p and e e the two components of the direction vector from the reference point to the calculated state. Note that these components remain constant as long as we move along the line connecting the reference and the trial state. Therefore, it is convenient to evaluate those components based on the trial state as With these quantities, the routine for the projected quantities on the yield is similar to a classical zero root finding algorithm for single variable functions, that is, finding 0 so that the resulting yield is identically 0: To that purpose, we make use of Brent's method (Brent (1971)) following the implementation presented in Press et al. (2007). The main advantage of our update procedure is that it permits to make use of different expressions for the yield and the flow laws.
As specified in section 2.4, the yield function can be undefined in certain regions of the stress invariants space (see, e.g., Stupkiewicz et al., 2014, which described this issue by means of the original Bigoni-Piccolroaz yield function; Bigoni & Piccolroaz, 2004). However, one can solve this issue by forcing the projections onto the yield surface as based on the square of the yield function in equation (34). For more details on such an update procedure, we invite the interested readers to the original paper by Stupkiewicz et al. (2014). The resulting projected values-calculated based on the square of the yield function-can be then used to evaluate the yield function given by equation (17) maintaining internal consistency of the units for the viscoplastic and damage flow rules.

Damage-Induced Weakening and Strain Localization
One of the main property controlling the onset of inelastic deformation and damage accumulation is the friction of the material. In this study, the yield function in the undamaged stress space is formulated in  Table 1 lists the physical properties used for this example.
terms of a constant friction coefficient. However, as already stated in the previous chapter, we can recast our damage formulation in terms of an effective friction coefficient as based on the effective stress tensor, the latter being a function of the damage intensity variable (Lyakhovsky et al., 1993(Lyakhovsky et al., , 2005Lyakhovsky & Ben-Zion, 2008) as In Figure 3 we illustrate the evolution of the effective friction coefficient as a function of the damage variable for different values of the elastic strain ratio. A general feature is an inverse state dependency of the friction coefficient on the level of accumulated damage. In addition, friction weakening is non-linear, with rates accelerating at higher values of the elastic strain ratio. As for Figure 3, the control of friction weakening/strengthening on the overall deformation mode is not constant but rather changes progressively during the evolution of the system, being intrinsically related not only to the actual damaged state, but also to the rate of damage accumulation experienced. This behavior resembles some of the basic features as derived from classical rate-and-state frictional models. However, in the current formulation the weakening/hardening mechanism depends exclusively on the free energy formalism adopted to derive the basic state equations to describe the evolution of the physical system (equation (8)).
In order to better illustrate and quantify the impact of damage-induced weakening on the evolution of localized deformation for low-porosity rocks we relied on a classical benchmark case for geodynamic modeling applications. We consider a setup as introduced by Kaus (2010), where a homogeneous lithosphere block is subjected to uniaxial compression. The main difference with respect to the classical setup is the rheology considered here, that is, in our study the material is governed by a damage visco-elasto-viscoplastic rheology. The initial friction of the material is 30 • and the inclusion at the bottom center of the domain is described as a Von Mises material. The geometry and boundary conditions adopted are illustrated in Figure 4. Upon compression, strain starts to localize around the weak inclusion leading to the formation of shear bands  which propagate to the top of the model. To better quantify the effects of damage rheology, we discuss two end-member cases in terms of damage viscosities . In the first realization, we impose a relative high damage viscosity, therefore forcing low damage accumulation in the system. In contrast, Case 2 has a lower damage viscosity, thereby promoting faster and higher magnitudes accumulation of damage. All properties used for this example are listed in Table 1.
In Figure 5 we plot the distribution of the second invariant of the total strain rate tensor (upper panels) as well as the distribution of the damage intensity variable (lower panels) after 0.5% of compression. For the case where damage accumulation is inhibited (Case 1), the effective friction coefficient does not undergo any sensitive variations in time. Therefore, strain localizations are controlled by the undamaged viscoplastic model. This results in a setup where maxima in strain rates coincide spatially with regions where the yield function reaches its maximum values, that is, along the edges of the developing bands, leading to a widening of the localized bands ( Figure 5a). The situation is different for Case 2. Higher levels of damage accumulation, significantly weakens the effective friction of the material. The main effect of the resolved weakening is to help maintaining the localized deformation within the initial shape of the bands, which are characterized by a higher degree of localization when compared to Case 1.

Pressure-Controlled Shear Bands Formation in Dilation and Compression Regimes
As already discussed in some detail in Jacquey and Cacace (2020), the origin and formation of shear bands is likely to be related to the presence of material heterogeneities and/or weak domains. Shear bands evolution and orientation, on the other hand, depend mainly on the (effective) confining pressure, which can be described in terms of the dilation angle. In the following application, we investigate the evolution of shear bands under pure shear by making use of a setup similar to the one introduced in Duretz et al. (2014Duretz et al. ( , 2018, but by systematically investigating the response of the system to different confining pressure. This is achieved by relying on a pressure-dependent capped yield, and by considering two end-member cases: (i) background pressure typical for dilation mode and (ii) background pressure in compression mode.
The geometry of the domain, set of boundary conditions as well as imposed confining pressure values for the two cases are illustrated in Figure 6. In both realizations, the material is considered to follow a visco-elasto-viscoplastic rheology with a capped yield as described in section 3.2. All properties are summarized in Table 2. Figure 7 illustrates the distribution of the inelastic porosity (related to the inelastic volumetric deformation) after (a) 1% of dilatant shear deformation and (b) 2% of compactant shear deformation. Dilatant behavior is characterized by an increase in the inelastic porosity, which localized in shear bands oriented at an angle of 33 • . Deformation is localized mainly within the evolving bands. The results obtained in compactant mode differs. In this case, deformation is more diffused within the whole domain area, though still forming regions of preferred localization in the forms of bands of decreased inelastic porosity. Those bands starts forming from the inclusion at the center of the model, and later evolves being oriented at 45 • around the inclusion.

Behavior Under Laboratory Conditions
Laboratory studies of porous rocks have demonstrated that the rock behavior in the inelastic and failure mode have some common attributes if investigated over a large range of loading conditions. In general,  a porous rock would deform inelastically by dilatancy under low confinement. Under these conditions, and, upon reaching a peak stress level, the sample would eventually fail in a macroscopic manner resulting in the development of shear bands and measurable strain softening. In contrast, under relatively high level of confinement, the same sample would preferentially undergo a decrease in its volume. Compactant rock behavior can be observed independently on the presence of any deviatoric stress. Indeed, compaction in porous rocks has been reported also under hydrostatic loading conditions, as related to an increase in the effective confining pressure. However, the volumetric response of a porous rocks under the same level of effective mean stress differs depending on the type of loading applied. Compaction is usually larger in samples undergoing non hydrostatic loading, a phenomenon known as shear-enhanced compaction (Baud et al., 2000;Curran & Carroll, 1979;Menéndez et al., 1996;Wong et al., 1992;Zhu et al., 1997). Dilatancy is usually considered as promoting strain localization (often predating the onset of macroscopic shear banding). However, deformation can also localized in compacting rocks.
Despite that the phenomenological behaviors of most porous rocks can be considered to be similar (localized brittle versus delocalized ductile), their bulk deformational response to pressure, temperature, strain rate, and presence of fluid can differ. In this section, we present a suite of numerical experiments under laboratory conditions. The setup is considered typical for triaxial mechanical loading of fully saturated samples under drained conditions. The test cases are intended to illustrate the influences of (i) the axial loading rate on the strength evolution, which controls the deformation history of the material; (ii) the effects of the (effective) confining pressure, and (iii) ambient temperature as they both give rise to the onset of semibrittle mechanical behavior. The aim of these applications is to showcase the ability of the damage poro-elasto-viscoplastic rheology to capture the basic response of porous rocks to varying forcing conditions as constrained by laboratory results. No attempt is here made to reproduce the details of a specific experiment, which would require a proper tuning of the basic rheological properties against available data and is subject of ongoing activities. Figure 8 sketches the basic ingredients of the conceptual model implemented to capture the transitional brittle-ductile behavior by means of the above described damage rheology. Given our overstress implementation, accounting for state variables dependency in the mechanical evolution requires to define a proper constitutive functional behavior of both the damage and viscoplastic viscosities. Such constitutive laws should enable to capture the basic characteristics of the deformation modes of porous rocks in the two end-member states. In the pure brittle field, deformation is followed by damage accumulation and shows a predominant dilatant behavior, being sensitive to the (effective) pressure. This experimental evidence is introduced in our model formulation by considering the damage viscosity to be an exponential function of the confining pressure leading to a decrease in damage rate upon an increase in confining pressure as ). In contrast, to a first order, ductile inelastic behavior can be considered not to be directly affected by variations in the pressure conditions, being only sensitive to temperature variations. It follows, that to reconcile this type of temperature sensitive behavior, the viscoplastic viscosity is described in the form of a classical Arrhenius's flow law ( p = 0 p exp ).
Laboratory experiments have demonstrated that in both end-member cases, deformation of porous rocks is highly sensitive to the rate of applied load. This aspect is implicitly taken into account in the rate-dependent (overstress) damage plastic model adopted, where the bulk mechanical response is a function of the resolved strain rate.
Simulations of triaxial loading of rock samples are performed on a three-dimensional mesh, including ∼20,000 hexagonal elements representing the cylindrical sample with a radius of 5 mm and a length of 20 mm. Following what stated above, we discuss different simulation runs where we systematically vary (i) the axial loading rate, (ii) the confining pressure, and (iii) the in situ temperature. For all simulations, we maintain the fluid pressure 10 MPa to represent drained conditions. Given the timescale of triaxial laboratory experiments, we neglect any viscoelastic effects in the following. Therefore, our numerical sample undergoes a poro-elasto-viscoplastic mechanical response to imposed loading. Table 3 lists all physical properties used. This set of property are considered representative of sandstone rock samples (Baud & Meredith, 1997;Brantut et al., 2013;Heap et al., 2009). Figure 9 illustrates the simulation results of the different triaxial tests by means of the evolution of the deviatoric stress (upper panels) and porosity change (lower panels) as a function of the axial strain.
The left graphs show the simulation results for different axial strain rate conditions and constant temperature and confining pressure. For all cases tested, the rock strength increases significantly with strain rate. higher magnitudes of strain rate (from 0.5 × 10 −5 to 1.5 × 10 −5 s −1 ) results in a net increase in the rock strength (approximately 100 MPa). At the lowest strain rate, the stress versus strain curve shows essentially a linear elastic behavior until the peak load. By increasing the strain rate, inelastic behavior is observed before the peak stress is reached at higher axial strain levels. Upon reaching the peak stress, all simulations undergo strain weakening (drop in rock strength), the latter being higher for higher strain rates. All three simulations are characterized by a porosity decrease during the initial stiffening phase. The amount of porosity reduction increases with increasing strain rate, thus implying a more efficient initial compaction at faster strain rates. At a later stage of loading, all cases show onset of dilatant behavior. Dilatancy starts before the rock reaches its peak strength, and accelerates during the strain weakening phase.
The middle panel shows the results obtained by simulating the response of the porous sample to different confining pressures at constant axial strain rate and temperature. The confining pressure was raised in three steps from 30 to 50 MPa and then to 90 MPa, corresponding to effective pressure loading of 20, 40, and 80 MPa, respectively. The maximum value of imposed confining pressure was intentionally kept below 100 MPa, in order to maintain the sample within the transitional brittle-ductile deformation behavior (as shown for the 90 MPa case). In all the three cases, the sample initially deforms elastically before reaching a critical stress value. Magnitudes of the critical stress increase with increasing confining pressure. Further compression leads to inelastic deformation, the onset of which also depends on the confining pressure. At low to moderate confining pressures (below 90 MPa), the stress-strain curves show a defined peak strength and a strain weakening postfailure phase until their residual strength. For the highest confining pressure, the rock does not show a distinctive peak strength and the work weakening phase is also smoother, indicative of a more ductile deformation behavior. This difference is also reflected in the porosity versus strain curves. The overall behavior for the sample subjected to confining pressure less than 90 MPa is similar to the one described above. Work hardening accompanied by compactant behavior, is followed by dilatancy before the peak strength and during the strain weakening phase. Dilation of the rock is significant at low confining pressures and leads to a considerable permanent volume increase (less than 1%) of the rock compared to their initial volume. Most of this volume increase occurs upon and after reaching the rock peak strength, that is in the strain weakening phase. In contrast, at the highest confining pressure, dilation substantially decreases and the porosity versus strain curve shows a more gradual compaction-to-dilation transition, an aspect that confirms the transition to a more ductile deformation mode at this level of confining pressure.
Increasing the in situ temperature from 25 up to 500 • C (right panel) leads to a a systematic reduction in strength. Depending on the temperature conditions, initial elastic behavior is followed by a different deformation mode. At low temperatures (25 • C) the rock gradually yields and hardens until its peak stress. The peak strength is then followed by strain weakening up to strains of about 1.5%. Strain hardening is characterized by compaction at small strain levels (less than 0.75%), gradually followed by a dilatant mode at larger strains. The onset of dilation is observed prior to the peak deviatoric stress. At intermediate temperatures (250 • C), the rock initially behaves elastically until it yields at around 150 MPa. Yielding leads to a smooth transition toward a monothonic weakening phase, which continues up to strains of about 2%, without residual strength being reached. Volume changes are indicative for compaction at strain up to 1.5%, followed by relative small amount of dilatation at larger strain. At the highest temperature the rock sample behaves in a more ductile fashion, reaching a constant deviatoric stress and a linear tendency for compaction at all strains.

Discussions
In this contribution we have described a physical framework and its numerical implementation to describe the mechanical response of porous rocks to varying loading conditions spanning the transition between brittle and ductile regimes. Brittle behavior is accounted for by a frictional damage weakening formulation, which depends both on the rate of loading and the state of the deformation. The coupling between yielding, damage of the material and its porosity evolution enables to account for both dilatant and compactant responses under differential loading. Damage accumulation and confining pressure exert a significant impact on the style of strain localization.
Based on the implemented formulations we can also reconcile typical material behaviors as observed during load constrained laboratory conditions. From the above results, we can conclude that the physical framework as described in this contribution allows to account for a range of material responses spanning from purely brittle at high strain rates, low confining pressure and low temperature to ductile increasing confining pressure and/or temperature. In addition, we also demonstrated how the evolution of the porosity during triaxial loading can be further used as an indicator for precursor activities to brittle failure (onset of dilation) or of enhanced ductile compaction. The add value is the scale independence of the theoretical and numerical framework, which enables a direct upscale from the laboratory to the field applications.
Rock deformation experiments remain the main source of observations by which to constrain the mechanical response of rock materials. Therefore, we have illustrated how our damage poro-elasto-viscoplastic rheology can be used to reproduce observations obtained during drained triaxial laboratory experiments at different conditions. One specific feature of interest, which has been extensively reported by recent experimental studies is the onset of porosity increase before reaching a peak in deviatoric stress. Dilatancy has also been demonstrated to lead to a power law increase of microfractures creation which, upon coalescence ultimately might result in macroscopic brittle failure of the whole sample (Arcangelis et al., 2016;Jones & Molnar, 1979;Renard et al., 2017;Schubnel et al., 2007). The onset of porosity increase can therefore be taken as an indication of substantial opening of microcracks compared to the overall elastic compression of the sample, which can in turn promote macroscopic shear failure along a major plane. The modeling results can capture this behavior, thought to be the main driving process of accelerating degradation of the material. Indeed, damage accumulation promotes dilation of the material which in term induces the two-stage weakening of the material: (i) the first being a weakening in terms of friction as illustrated in section 4.1 and (ii) the second related to an overall change in the shape of the capped yield function. The critical pressure capping the yield depends on the actual porosity. Therefore, an increase in porosity will induce a decrease in the critical pressure, thus accelerating damage accumulation and inelastic strain. This self-promoting weakening processes-initiated at the onset of dilatancy-results in a power law increase of damage in a similar fashion to Omori's law often used in earthquake's physics and provide a quantitative run away point for instability to occur, as also observed during laboratory experiments.
We have investigated in this contribution the impacts of a damage poro-visco-elasto-viscoplastic formulation giving rise to semibrittle, semiductile deformation mechanisms on the initiation and evolution of strain localization at different timescales. While the outcomes of this study give insights into how localized features such as fractures and faults may generate in porous rocks, alternative formulations may be more suited to study the mechanical behavior of existing faulted structures. Those include models describing semibrittle deformation relying on a rate-and-state friction formulation linking the shear and normal stresses to the fault plane (Barbot & Fialko, 2010;Barbot, 2019;Goswami & Barbot, 2018;Noda & Shimamoto, 2010;Shimamoto & Noda, 2014). While these alternative formulations may differ in some aspect to the more continuum-based study presented here, Lyakhovsky et al. (2005) demonstrated that the damage rheology adopted in this study (based on the same Helmholtz free energy) exhibits some properties of a rate-and-state friction formulation.
Future work will focus on a proper upscaling of the formulation to natural and engineering (reservoir) conditions. Loading rates that can be reached in the laboratory are difficult to reconcile with those found in either natural tectonic settings and/or engineered reservoirs. Typical axial loading rates during triaxial experiments are in the order of 1.0 × 10 −5 s −1 whereas tectonic strain rates are in the order of 1.0 × 10 −15 s −1 . Such an upscaling can be achieved by considering the strain rate dependency of the rheological framework proposed in this study.
In addition, creep experiments can be also used to get further insights in the precursor dynamics, especially when applied to low strain, intracontinental regions. If supported by (still scarce) field observations, such creep tests can provide proper bounds to the material properties and their evolution and their impact on the mechanical response of these natural systems. Integrating laboratory observations of triaxial experiments at high temperatures (Violay et al., 2012(Violay et al., , 2017 and under creep conditions (Brantut et al., 2013;Heap et al., 2009Heap et al., , 2011 is the subject of ongoing work.

Conclusion
In this contribution, we have introduced a novel physical framework for modeling thermo-hydromechanical coupled processes driving the semi-brittle, semi-ductile behavior of porous rocks. We make use of a damage rheology to introduce material weakening driving brittle failure, which is then coupled to a thermodynamically consistent formulation of the porous response upon inelastic deformation. We have illustrated the impact of a damage poro-visco-elasto-viscoplastic rheology on the style of localized deformation in different pressure conditions and demonstrated the reliability of the predicted material responses at the laboratory scale under varying strain rate, confining pressure, and temperature conditions. Improving our understanding of semibrittle, semibrittle deformation for proper targeting of high-enthalpy unconventional reservoirs or mitigation of geological hazards-either induced or natural-requires properly identifying the driving forces and stability conditions inherent to the individual deformation mechanisms. The main conclusions that we derive can be summarized in the following points: • Material weakening due to degradation of the pore grain structure is an essential though usually neglected component to fully capture brittle deformation in numerical experiments. • Porosity evolution plays a major role in controlling the mechanical response of porous rocks upon loading likely providing indications for precursor activities leading to macroscopic material failure. • Despite ongoing improvements (as pointed out in the final discussion), the physical framework described in the two companion papers can already help to advance the current understanding of semi-brittle, semi-ductile deformation occurring across different scales, from the laboratory to the field.