A new anisotropic poroelasticity model to describe damage accumulation during cyclic triaxial loading of rock

SUMMARY

Crustal rocks undergo repeated cycles of stress over time. In complex tectonic environments where stresses may evolve both spatially and temporally, such as volcanoes or active fault zones, these rocks may experience not only cyclic loading and unloading, but also rotation and/or reorientation of stresses. In such situations, any resulting crack distributions form sequentially and may therefore be highly anisotropic. Thus, the tectonic history of the crust as recorded in deformed rocks may include evidence for complex stress paths, encompassing different magnitudes and orientations. Despite this, the ways in which variations in principal stresses influence the evolution of anisotropic crack distributions remain poorly constrained. In this work, we build on the previous non-linear anisotropic damage rheology model by presenting a newly developed poroelastic rheological model which accounts for both coupled anisotropic damage and porosity evolution. The new model shares the main features of previously developed anisotropic damage and scalar poroelastic damage models, including the ability to simulate the entire yield curve through a single formulation. In the new model, the yield condition is defined in terms of invariants of the strain tensor, and so the new formulation operates with directional yield conditions (different values for each principal direction) depending on the damage tensor and triaxial loading conditions. This allows us to discern evolving yield conditions for each principal stress direction and fit the measured amounts of accumulated damage from previous loading cycles. Coupling between anisotropic damage and anisotropic compaction along with the damage-dependent yield condition produces a reasonable fit to the experimentally obtained stress–strain curves. Furthermore, the simulated time-dependent cumulative damage is well correlated with experimentally observed acoustic emissions during cyclic loading in different directions. As such, we are able to recreate many of the features of the experimentally observed directional 3-D Kaiser ‘damage memory’ effect.

Creep and deformation, Elasticity and anelasticity, Fracture and flow, Mechanics, theory, and modelling

1 INTRODUCTION

It is well-established that crack damage is generated in brittle rocks that are subjected to a level of stress above some crack initiation threshold, and that this cracking results in the output of elastic wave energy in the form of acoustic emissions (AEs) (e.g. Meredith et al. 1990; Holcomb 1993; Lockner 1993a). During cyclic loading, cracks close elastically during unloading and re-open elastically during reloading. If the level of stress during reloading remains below the peak stress level attained in any previous loading cycle, then no new cracking occurs and no further cracking-related AE is generated. However, on any loading cycle where the previous peak stress is reached or exceed, new cracks are formed and are accompanied by concomitant AE output (Kurita & Fujii 1979; Holcomb & Costin 1986; Li & Nordlund 1993; Lockner 1993a; Pestman & Van Munster 1996; Lavrov 2001, 2003; Browning et al. 2017, 2018). This observation of AE output only recommencing when the previous maximum stress level is exceeded is known as the Kaiser effect (Kaiser 1953) and is related to the ability of a material to accumulate and reproduce information about previously experienced stress states. However, most experiments that have probed aspects of the Kaiser effect to date have been conducted during either uniaxial or conventional triaxial compression experiments and so have not been able to probe fully for any directionality in crack damage accumulation related to the orientation of principal stresses.

More recently, Browning et al. (2017, 2018) investigated the occurrence of a Kaiser effect in samples of Darley Dale sandstone subjected to both conventional and true triaxial stress conditions. Samples were loaded sequentially to increasing levels of peak stress, both with the maximum principal stress maintained in the same orientation and with the maximum principal stress rotated and applied sequentially in three orthogonal orientations. Their results showed that, under true triaxial loading, crack damage is a distinctly directional phenomenon, such that rocks can exhibit a 3-D, directionally dependent Kaiser effect, with AE only being generated when the previous peak stress in any specific orientation was exceeded. They therefore concluded that the Kaiser effect should more accurately be described as a damage memory effect rather than a stress memory effect.

Traditionally, the analysis of rock deformation and failure criteria has been formulated by, for example, a classical Coulomb–Mohr condition that defines brittle failure, and by a yield cap criterion that defines cataclastic flow (e.g. Issen & Rudnicki 2000). However, these formulations usually ignore any connection between yield stress and the amount of inelastic damage in the form of microcracks, voids or other flaws that leave the yield stress unchanged, and in doing so, ignore the underlying principle of the Kaiser effect. Laboratory experiments on porous rocks demonstrate evidence of overall strain hardening and yield cap growth attributed to plasticity and porosity loss (Baud et al. 2006; Tembe et al. 2008; Bedford et al. 2018). Several models have been developed for elasto-plastic deformation of isotropic soils, which are commonly formulated in a framework of continuum mechanics and can be successfully applied to model rock behaviour with complex yield conditions.

For example, the original Cam-Clay model (Roscoe & Burland 1968) provides a description for the stress versus inelastic strain behaviour for yield envelopes of any shape defined in stress space (Muir Wood 1990). Modified Cam‐Clay yield functions were successively used in geo-mechanical modelling of hydrocarbon reservoirs (Chan et al. 2004; Crawford et al. 2011) and in more generic studies of inelastic sandstone deformation (Schultz & Siddharthan 2005; Skurtveit et al. 2013). In the modified Cam‐clay formulation, the stress conditions required for yield are described by the elliptical function of differential and mean effective stress values. Grueschow & Rudnicki (2005) discussed the various models that incorporate different shapes of the evolving yield caps and compared their model with previous suggested models by DiMaggio & Sandler (1971) and Carroll (1991). These studies demonstrated that inelastic behaviour of porous rocks are well described by various plasticity models. Pijnenburg et al. (2019) quantified the elastic and inelastic contributions to the total deformation behaviour of Slochteren sandstones and concluded that not only the expanding yield envelopes, but also change in the elastic moduli should be considered in order to obtain a proper fit to the experimental stress–strain data. Damage rheology models are able to incorporate changes in both the local elastic properties and the form of the porosity-induced yield cap such that deformation patterns and modes of failure can be analysed alongside the yield cap growth (e.g. Bercovici et al. 2001; Stefanov et al. 2011; Lyakhovsky et al. 2015; Vorobiev 2019).

It has been suggested that the observed Kaiser effect in rocks indicates strain hardening in consecutive cycles such that the phenomena can be attributed to changes in yield surface due to damage accumulation (Holcomb 1993). Damage and yield surface growth are then likely coupled, and constraints on this coupling can aid interpretation of damage localization patterns and the Kaiser effect (Gajst et al. 2020). Damage evolution and time-dependent behaviour in low porosity sandstones have been investigated (Choens et al. 2021) through creep and conventional triaxial experiments, and numerical analyses. Such quasi-static and creep experiments have been successfully simulated using the modified poroelastic damage model of Lyakhovsky et al. (2015). As damage accumulates in the samples, the yield cap evolves to keep pace with the strain accumulation.

True triaxial experiments reported by Browning et al. (2017, 2018) demonstrated that the orientation of distributed microcracks in Darley Dale sandstone samples are essentially anisotropic and therefore require an extension of isotropic damage models using a scalar damage parameter and a more complex formulation that introduces a second-order damage tensor (Panteleev et al. 2021). The goal of this paper is, therefore, to provide a complete quantitative description of the rheological model with directional yield conditions (i.e. different values for each principal direction) depending on the damage tensor and triaxial loading conditions. The new model combines and extends the results of the previously developed anisotropic damage model of Panteleev et al. (2021) and the scalar poroelastic damage model of Lyakhovsky et al. (2015). The new analysis includes the ability to simulate yield curves through a single formulation and recreates many of the features of the experimentally observed directionally dependent Kaiser damage memory effect reported by Browning et al. (2018).

2 ANISOTROPIC POROELASTIC DAMAGE MODEL

2.1 Damage and porosity

Rock deformation is associated with the formation and growth of internal flaws. From a mechanical point of view, these flaws can be divided into two classes: (1) microcracks (damage) contained in the matrix of a porous rock which act as primary stress raisers or stress concentrators and hence contribute to brittle failure and (2) pores which can deform and before their collapse act to dissipate or accommodate stress and hence contribute to distributed flow. For an isotropic rock with a sufficiently large number of microcracks and pores, one can define a representative volume in which the flaw density is uniform and described by two scalar variables, damage (⁠|${\rm{\alpha }}$|⁠) and porosity (⁠|${\rm{\varphi }}$|⁠). The damage variable is a mechanical variable, which is responsible for the change in material stiffness and brittle failure at a critical level of damage. For anisotropic rocks, we can consider a damage tensor, |${{\rm{\Omega }}_{ik}}$|⁠, which represents not only the density of microcracks, but also their orientations. The porosity variable is a geometrical property representing the volume fraction of pores during and after deformation. As an alternative to porosity, we define a compaction-strain tensor, |${{\rm{\psi }}_{ij}}$|⁠, which is equal to the accumulated irreversible strain resulting from loading and unloading. This tensor then represents not only the pore volume change, but also deviations in the shape of the pores.

In the following subsections we describe the general thermodynamic approach used to construct the scalar damage and poroelastic damage model, and provide the main equations of the new anisotropic poroelastic damage model. Detailed thermodynamic relations are provided in Appendix  A and specific relation for the isotropic and anisotropic model formulations are provided in Appendixes  B and  C.

2.2 General thermodynamic approach

We derive the main equations of the poroelastic damage model using the basic relations of irreversible thermodynamics, which provide constraints on the rates of dissipative processes (e.g. Onsager 1931; Biot 1955; Prigogine 1955; Truesdell & Noll 2004; DeGroot & Mazur 2013). This approach has been applied successfully to understand the kinetics of chemical reactions and phase transitions (e.g. Fitts 1962; DeGroot & Mazur 2013), and as the basis for variational methods of continuous media models (e.g. Sedov 1968; Malvern 1969, 1997; Berdichevsky 2009). The constitutive behaviour of the material, and flow rules controlling the kinetics of related irreversible processes, is then entirely defined by specification of two potentials. The first is the free energy, F, and the second is the dissipation function or local entropy production, Γ. This approach has been used as the basis for other damage models (e.g. Valanis 1990; Hansen & Schreyer 1994; Lyakhovsky et al. 1997; Bercovici et al. 2001; Hamiel et al.2004a; Gaede et al. 2013). Following Onsager (1931), who theoretically generalized the empirical laws of Fourier, Ohm, Fick and Navier (see review by Martyushev & Seleznev 2006), we represent the specific local entropy production as a product of thermodynamic fluxes and thermodynamic forces. For small deviations from equilibrium, the Onsager principle can be obtained from the maximum entropy production principle, the maximum dissipation rate of mechanical energy, or the von Mises principle (e.g. Martyushev & Seleznev 2006; Ziegler 2012). We now discuss the different forms of the energy function, beginning with the scalar isotropic damage formulation, then the coupling isotropic damage and porosity model, and finally we formulate the anisotropic model. The energy and entropy balance equations and general thermodynamic relations are provided in Appendix  A.

2.3 Scalar damage and poroelastic damage model

The free energy of a solid (F) in the local damage model of Lyakhovsky et al. (1997) is assumed to be a function of the state variables, which are the temperature T, the elastic strain tensor |${\varepsilon _{ij}} = g_{ij}^{( t )}\ \ - g_{ij}^{( 0 )}$| (the difference between the total strain tensor |$g_{ij}^{( t )}$| and the irreversible strain tensor |$g_{ij}^{( 0 )}$|⁠), and the scalar damage variable α:

$$\begin{eqnarray*} F\ = \ F\left( {T,{{\rm{\varepsilon }}_{ij}},{\rm{\alpha }}} \right). \end{eqnarray*}$$

(1)

Using the balance equations for the energy and entropy, the Gibbs relation and the Murnaghan (1937) definition of the stress tensor, part of local entropy (⁠|$\Gamma $|⁠) production associated with evolving damage is:

$$\begin{eqnarray*} \Gamma \ = \ - \frac{{\partial F}}{{\partial {\rm{\alpha }}}}\frac{{{\rm d}{\rm{\alpha }}}}{{{\rm d}t}} \ge 0. \end{eqnarray*}$$

(2)

The complete thermodynamic derivations are presented in Appendix  A, where all the dissipation processes are fully discussed. Following the Onsager (1931) principle, the kinetic relation for damage evolution is:

$$\begin{eqnarray*} \frac{{{\rm d}{\rm{\alpha }}}}{{{\rm d}t}} = \ - C\ \frac{{\partial F}}{{\partial {\rm{\alpha }}}}, \end{eqnarray*}$$

(3)

where C is the positive kinetic coefficient, which may be either constant or depend on the state variables.

Hamiel et al. (2004b) and then Lyakhovsky et al. (2015) extended the scalar damage model to permit coupling of damage and porosity in their formulations. They followed Biot's theory of poroelasticity (Biot 1941, 1956) representing the free energy of a poroelastic medium as a sum of the elastic energy, and the Biot poroelastic coupling terms of the saturated medium. The free energy (1) is extended to be a function of two additional state variables, fluid volume content, ζ, and material porosity, φ:

$$\begin{eqnarray*} F\ = \ F\left( {T,{{\rm{\varepsilon }}_{ij}},{\rm{\alpha }},{\rm{\varphi }},{\rm{\zeta }}} \right). \end{eqnarray*}$$

(4)

As both porosity and damage can evolve with time during deformation, their coupled kinetic equations are derived using a similar balance equation leading to the following local entropy production (Hamiel et al. 2004b, Lyakhovsky et al. 2015):

$$\begin{eqnarray*} \Gamma \ = \ - \frac{{\partial F}}{{\partial {\rm{\alpha }}}}\frac{{{\rm d}{\rm{\alpha }}}}{{{\rm d}t}} - \left( {\frac{{\partial F}}{{\partial {\rm{\varphi }}}} + {\sigma _m}} \right)\frac{{{\rm d}{\rm{\varphi }}}}{{{\rm d}t}} \ge 0, \end{eqnarray*}$$

(5)

where |$\ {\sigma _m} = \ - {\sigma _{kk}}/3$| is the mean stress. Once more, adopting the relations from Onsager (1931) gives a set of two coupled differential equations (Malvern 1969; DeGroot & Mazur 2013) which define the damage and porosity evolution:

$$\begin{eqnarray*} \frac{{{\rm d}{\rm{\varphi }}}}{{{\rm d}t}} = \ - {C_{{\rm{\varphi \varphi }}}}\left( {\frac{{\partial F}}{{\partial {\rm{\varphi }}}} + {\sigma _m}} \right) - {C_{{\rm{\varphi \alpha }}}}\frac{{\partial F}}{{\partial {\rm{\alpha }}}} \end{eqnarray*}$$

(6a)

$$\begin{eqnarray*} \frac{{{\rm d}{\rm{\alpha }}}}{{{\rm d}t}} = \ - {C_{{\rm{\alpha \varphi }}}}\left( {\frac{{\partial F}}{{\partial {\rm{\varphi }}}} + {\sigma _m}} \right) - {C_{{\rm{\alpha \alpha }}}}\frac{{\partial F}}{{\partial {\rm{\alpha }}}}. \end{eqnarray*}$$

(6b)

These phenomenological kinetic equations guarantee the non-negative value of entropy production if the matrix of the kinetic coefficients;

$$\begin{eqnarray*} {C_{ij}} = \left| {\begin{array}{@{}*{2}{c}@{}} {{C_{{\rm{\varphi \varphi }}}}}&{{C_{{\rm{\varphi \alpha }}}}}\\ {{C_{{\rm{\alpha \varphi }}}}}&{{C_{{\rm{\alpha \alpha }}}}} \end{array}} \right| \end{eqnarray*}$$

(7)

meets the following conditions (Malvern 1969; DeGroot & Mazur 2013): the diagonal cells (⁠|${C_{{\rm{\varphi \varphi }}}}$|⁠, |${C_{{\rm{\alpha \alpha }}}}$|⁠) must be positive, and the off-diagonal terms are usually taken to be either symmetric or antisymmetric. Following the poroelastic damage model of Hamiel et al. (2004b) and Lyakhovsky et al. (2015) we adopt an antisymmetric structure (⁠|${C_{{\rm{\varphi \alpha }}}} = \ - \ \ {C_{{\rm{\alpha \varphi }}}} = \ D$|⁠) of the kinetic matrix (7). These conditions assure positive dissipation, as in eq. 5. Larger D-values then lead to an earlier onset of damage and enhanced accumulation under the same confinement conditions. Hamiel et al. (2004b) and Lyakhovsky et al. (2015) discussed slightly different forms of the stress- or strain-dependent D-value and demonstrated how their scalar poroelastic model reproduces a yield cap and its evolution (see Appendix  B for details). Recently Gajst et al. (2020) suggested a model with exponential damage-dependent D-value:

$$\begin{eqnarray*} D\left( \alpha \right)\ = {D_1}\ \ {e^{ - {D_2}\alpha }}\ {\left( { - {I_1}} \right)^N}\sqrt {{I_2}} , \end{eqnarray*}$$

(8)

where the first coefficient D₁ stands for the initial D-value of the damage-free material, and the second coefficient D₂ controls its decrease with increasing damage. The role of the exponent N > 1 is to control the shape of the yield cap; and this is further discussed in the supplementary materials. Gajst et al. (2020) demonstrated that the decrease in the D-value with accumulated damage shifts the yield condition or onset of damage to higher stress values and successfully reproduces the Kaiser effect. Since the model is formulated in terms of scalar damage and strain invariants, it accurately reproduces the isotropic Kaiser damage-memory effect, but does not consider the effect of microcrack orientation and stress rotation. The experimentally observed directionally dependent Kaiser damage-memory effect (Browning et al. 2018) hence requires an anisotropic formulation.

2.4 Anisotropic poroelastic damage model

Recently Panteleev et al. (2021) extended the scalar isotropic damage model by developing a theoretical model for materials with orthotropic symmetry which describes the material damage using a second rank symmetric tensor, |${{\rm{\Omega }}_{ik}}$|⁠, in which the principal directions match the orientation of the principal loading axes. This assumption is supported by results of true triaxial experiments (Browning et al. 2017, 2018) which demonstrated that the orientation of distributed microcracks was related to the level and orientation of the principal stresses. Therefore, most of the equations that follow are written with respect to the principal loading directions and stresses, while the complete 3-D formulation is presented in Appendix  C.

The scalar damage variable α in the free energy form of eq. (1) and the poroelastic model with the energy form from eq. (4) is substituted by a damage tensor, |${{\rm{\Omega }}_{ik}}$|⁠. For the case of an isotropic material ( |${{\rm{\Omega }}_{ij}} = \ {\rm{\Omega }}\ {\delta _{ij}}$|⁠), the anisotropic formulation reduces to the scalar model with the damage |$\alpha $| equal to a squared value, |$\alpha \ = \ {{\rm{\Omega }}^2}$|⁠. In addition, Lyakhovsky et al. (2022) showed that the deformation of pore space is inherently 3-D and, as such, the compaction-strain strain tensor, |${{\rm{\psi }}_{ij}}$|⁠, should replace porosity in the governing equations. The suggested energy function includes these two tensor state variables; the damage tensor |${{\rm{\Omega }}_{ik}}$|⁠, and the compaction-strain tensor |${{\rm{\psi }}_{ij}}$|⁠:

$$\begin{eqnarray*} F\ = \ F\left( {T,{{\rm{\varepsilon }}_{ij}},{{\rm{\Omega }}_{ij}},{{\rm{\psi }}_{ij}},{\rm{\zeta }}} \right). \end{eqnarray*}$$

(9)

The elastic strain tensor |${\varepsilon _{ij}} = g_{ij}^{( t )}\ \ - {{\rm{\psi }}_{ij}}$| is now defined as the difference between the total strain tensor |$g_{ij}^{( t )}$| and the tensor |${{\rm{\psi }}_{ij}}$|⁠. The diagonal part of this tensor, |$\varphi \ = \ {{\rm{\psi }}_{ij}}{\delta _{ij}}$| represents the material porosity, while the deviatoric components (⁠|${{\rm{\psi }}_{ij}} - \frac{1}{3}\varphi \ {\delta _{ij}}$|⁠) are associated with anisotropic compaction and other mechanisms related to the irreversible strain accumulation. Using the energy form (9), the dissipation associated with evolving tensors |${{\rm{\Omega }}_{ij}},{{\rm{\psi }}_{ij}}$|⁠, consists of two terms which are proportional to the their rate of change (see Appendix  C for detailed derivations):

$$\begin{eqnarray*} \Gamma \ = \ - \frac{{\partial F}}{{\partial {\Omega _{ij}}}}\frac{{{\rm d}{\Omega _{ij}}}}{{{\rm d}t}} + \left( {{{\rm{\sigma }}_{ij}} - \frac{{\partial F}}{{\partial {\psi _{ij}}}}} \right)\frac{{{\rm d}{\psi _{ij}}}}{{{\rm d}t}} > 0. \end{eqnarray*}$$

(10)

The phenomenological kinetic equations share the same structure with the poroelastic formulation (6), but connect the tensor quantities as follows:

$$\begin{eqnarray*} \frac{{{\rm d}{\psi _{ij}}}}{{{\rm d}t}} = C_{ijnm}^{\psi \psi }\ \ \left( {{{\rm{\sigma }}_{ij}} - \frac{{\partial F}}{{\partial {\psi _{ij}}}}} \right) - C_{ijnm}^{\psi \Omega }\frac{{\partial F}}{{\partial {\Omega _{ij}}}} \end{eqnarray*}$$

(11a)

$$\begin{eqnarray*} \frac{{{\rm d}{\Omega _{nm}}}}{{{\rm d}t}} = C_{ijnm}^{\Omega \psi }\ \ \left( {{{\rm{\sigma }}_{ij}} - \frac{{\partial F}}{{\partial {\psi _{ij}}}}} \right) - C_{ijnm}^{\Omega \Omega }\frac{{\partial F}}{{\partial {\Omega _{ij}}}}. \end{eqnarray*}$$

(11b)

Instead of the matrix (7) of the kinetic coefficients, every matrix term becomes a fourth-rank tensor that can be written as:

$$\begin{eqnarray*} {C_{ijnm}} = \left| {\begin{array}{@{}*{2}{c}@{}} {C_{ijnm}^{\psi \psi }}&{C_{ijnm}^{\psi \Omega }}\\ {C_{ijnm}^{\Omega \psi }}&{C_{ijnm}^{\Omega \Omega }} \end{array}} \right|{\rm{\ }}. \end{eqnarray*}$$

(12)

The kinetic eq. (11) guarantee a non-negative value of entropy production if the cells of the matrix of the kinetic coefficients meet conditions like those of the poroelastic model: (1) matrices of the diagonal cells (⁠|$C_{ijnm}^{\psi \psi }$|⁠, |$C_{ijnm}^{\Omega \Omega }$|⁠) must be positively defined and (2) we also adopt an antisymmetric structure ( |$C_{ijnm}^{\Omega \psi } = \ - C_{ijnm}^{\psi \Omega }$|⁠) for the off-diagonal terms, as was done previously for the poroelastic model.

In the next section, we specify the energy function (9) and kinetic coefficients (12), and then demonstrate the main model features.

2.5 Energy function and kinetic equations, anisotropic model

The energy function for the anisotropic damage model includes a damage tensor |${\Omega _{ij}}$| and so cannot be formulated only in terms of invariants of the strain tensor:

$$\begin{eqnarray*} \begin{array}{@{}*{1}{c}@{}} {\ \ {I_1} = {\rm{\ }}{{\rm{\varepsilon }}_{ij}}{{\rm{\delta }}_{ij}}}\\ {\ \ {I_2} = {\rm{\ }}{{\rm{\varepsilon }}_{ij}}{{\rm{\varepsilon }}_{ij}}} .\end{array} \end{eqnarray*}$$

(13)

Following Murti et al. (1991), Zhang & Cai (2010) and Panteleev et al. (2021) incorporated invariants |$I_1^{( {\rm{\Omega }} )}$| and |$I_2^{( {\rm{\Omega }} )}$| of the tensor |$\varepsilon _{ij}^{( {\rm{\Omega }} )} = \frac{1}{2}\ \ ( {{\varepsilon _{ik}}{{\rm{\Omega }}_{kj}} + {\varepsilon _{jk}}{{\rm{\Omega }}_{ki}}} )$|⁠. In the coordinate system of the principal damage values, these invariants are (see Appendix  C for the general case):

$$\begin{eqnarray*} &&{I_1^{\left( {\rm{\Omega }} \right)} = {\varepsilon _{11}}{{\rm{\Omega }}_1} + {\varepsilon _{22}}{{\rm{\Omega }}_2} + {\varepsilon _{33}}{{\rm{\Omega }}_3}}\nonumber \\ && {I_2^{\left( {\rm{\Omega }} \right)} = {{\left( {{\varepsilon _{11}}{{\rm{\Omega }}_1}} \right)}^2} + {{\left( {{\varepsilon _{22}}{{\rm{\Omega }}_2}} \right)}^2} + {{\left( {{\varepsilon _{33}}{{\rm{\Omega }}_3}} \right)}^2} + }\nonumber \\ && {\frac{1}{2}\varepsilon _{12}^2{{\left( {{{\rm{\Omega }}_1} + {{\rm{\Omega }}_2}} \right)}^2} + \frac{1}{2}\varepsilon _{13}^2{{\left( {{{\rm{\Omega }}_1} + {{\rm{\Omega }}_3}} \right)}^2} + \frac{1}{2}\varepsilon _{23}^2{{\left( {{{\rm{\Omega }}_2} + {{\rm{\Omega }}_3}} \right)}^2}}. \end{eqnarray*}$$

(14)

We extend the energy function of Panteleev et al. (2021) using additional terms of Biot's theory of poroelasticity (Biot 1941, 1956; see also Hamiel et al. 2004b and Lyakhovsky et al. 2015):

$$\begin{eqnarray*} F\ = \left[ {\frac{{{\lambda _0}}}{2}I_1^2 + {\mu _0}{I_2}} \right]\ \ + \left[ {{\mu _1}I_2^{\left( {\rm{\Omega }} \right)} - \gamma I_1^{\left( {\rm{\Omega }} \right)}\sqrt {I_2^{\left( {\rm{\Omega }} \right)}} } \right] + \frac{M}{2}{\left[ {\beta {I_1} - \zeta + {{\rm{\psi }}_{ij}}{\delta _{ij}}} \right]^2} + Ch\frac{\gamma }{2}{{\rm{\Omega }}_{ij}}{{\rm{\Omega }}_{ij}}. \end{eqnarray*}$$

(15)

The energy function for non-linear poroelastic damaged media includes two Hookean terms with the Lamé drained moduli of the intact (damage-free) rock |${\lambda _0},\ {\mu _0}$| and two second order terms with strain invariants |$I_1^{( {\rm{\Omega }} )}$| and |$I_2^{( {\rm{\Omega }} )}$|⁠. The modulus |${\mu _1}$| controls the reduction of the effective shear modulus, and the coupling modulus γ is responsible for enhanced non-linearity with damage accumulation (Panteleev et al. 2021). The third term in squared brackets, with Biot modulus M and coefficient |$\beta ,$| differs from the classical poroelasticity only by the term |${{\rm{\psi }}_{ij}}{\delta _{ij}},$| that represents the porosity. Similarly to the scalar model (Appendix  B), we introduce the damage-dependent term with the coefficient Ch (Gajst 2020), which allows us to account for the cohesive forces that influence rock fracture under low confining pressures.

Following the definitions of the scalar damage model, we use |${\mu _1} = {\xi _0}\ \ \gamma $| with critical ratio |${\xi _0} = {I_1}\ /\sqrt {{I_2}} $| of the strain invariants (13). The |${\xi _0}$| value is related to the internal friction angle of the intact rock (Agnon & Lyakhovsky 1995) and controls the onset of damage accumulation in the scalar damage model as well as in the anisotropic model for the material with isotropic damage ( |${{\rm{\Omega }}_{ij}} = \ {\rm{\Omega }}\ {\delta _{ij}}$|⁠). With this notation, the stress–strain constitutive relation for for k-component of the principal stress and damage values is (see Appendix  C for complete relation):

$$\begin{eqnarray*} {\sigma _k} = {\lambda _0}\ \ {I_1} + 2{\mu _0}{\varepsilon _k} - \gamma {{\rm{\Omega }}_k}\sqrt {I_2^{\left( {\rm{\Omega }} \right)}} + \gamma \left( {2{\xi _0} - \frac{{I_1^{\left( {\rm{\Omega }} \right)}}}{{\sqrt {I_2^{\left( {\rm{\Omega }} \right)}} }}} \right){\varepsilon _k}{\rm{\Omega }}_k^2 + \beta M{\rm{\ }}\left( {\beta {I_1} - \zeta + {{\rm{\psi }}_k}} \right) \end{eqnarray*}$$

(16)

and the fluid pressure is:

$$\begin{eqnarray*} {p_f} = \frac{{\partial F}}{{\partial \zeta }}\ \ = \ M\ \left( { - \beta {I_1} + \zeta - {{\rm{\psi }}_n}} \right). \end{eqnarray*}$$

(17)

Similarly to the Biot poroelasticity, the effective stress is defined as:

$$\begin{eqnarray*} \ \ \sigma _{ij}^{\rm eff} = \ {\sigma _{ij}} - \frac{{\partial F}}{{\partial {{\rm{\psi }}_{ij}}}}. \end{eqnarray*}$$

(18)

Kinetic coefficients (12) should be defined in order to provide the complete form of the kinetic eq. (11). The matrix |$C_{ijnm}^{\Omega \Omega }$| multiplied by |$\partial F/\partial {\Omega _{ij}}$| defines the damage accumulation rate, driven by the thermodynamic force associated with the damage-dependent energy change. The most conservative assumption to define the components of this matrix is the absence of any interaction between different damage components on their kinetics:

$$\begin{eqnarray*} C_{ijnm}^{\Omega \Omega } = \frac{L}{{\sqrt {{{\rm{\Omega }}_{ij}}{{\rm{\Omega }}_{ij}}} }}\ \ \ \left( {{\delta _{ik}}\ {\delta _{jn}} + {\delta _{in}}\ {\delta _{jk}}} \right). \end{eqnarray*}$$

(19)

This form of damage kinetics was verified by Panteleev et al. (2021) using results from true triaxial rock mechanics experiments, and is therefore adopted here for the poroelastic model.

The off-diagonal antisymmetric coefficient, noted here as |${D_{ijkn}} = C_{ijnm}^{\Omega \psi }\ \ = \ - C_{ijnm}^{\psi \Omega },$| controls the coupling between irreversible strain (porosity) and damage accumulation. Extending the Gajst et al. (2020) model with an exponential damage-dependent D-value (8) to the tensor form and using the same type of strain dependency, we suggest the following form of the coupling kinetic coefficient |${D_{ijkn}}$|⁠:

$$\begin{eqnarray*} {D_{ijkn}} = {D_1}\ \ \exp\left[ { - {D_2}{{\rm{\Omega }}_{ij}}} \right]\ {\delta _{kn}}\ {\left( { - {I_1}} \right)^N}\sqrt {{I_2}} , \end{eqnarray*}$$

(20)

where we use the standard definition of the exponent of the tensor |${{\boldsymbol{X}}_{ij}}$| by means of its series representation (Hirsch et al. 1974):

$$\begin{eqnarray*} \exp\ \ \left( {{{\boldsymbol{X}}_{ij}}} \right) = \ \mathop \sum \limits_{m\ = \ 0}^\infty \frac{{{\boldsymbol{X}}_{ij}^m}}{{m!}}. \end{eqnarray*}$$

(21)

Note that the principal values of the tensor |$\exp( {{{\boldsymbol{X}}_{ij}}} )$| are equal to the exponent of the principal value |$\exp( {{{\boldsymbol{X}}_k}} )$|⁠.

Given the kinetic coefficients (19), (20), the equation for the damage evolution (11b) for principal components becomes (see Appendix  C for complete relation):

$$\begin{eqnarray*} \frac{{{\rm d}{\Omega _k}}}{{{\rm d}t}} &=& K\ {D_1}\exp\left[ { - {D_2}{{\rm{\Omega }}_k}} \right]\ {\left( { - {I_1}} \right)^{N + 1}}\sqrt {{I_2}} \ + \nonumber \\&& \frac{L}{{\sqrt {{{\rm{\Omega }}_l}{{\rm{\Omega }}_l}} }}\left[ {{\varepsilon _k}\sqrt {I_2^{\left( {\rm{\Omega }} \right)}} + \left( {\frac{{I_1^{\left( {\rm{\Omega }} \right)}}}{{\sqrt {I_2^{\left( {\rm{\Omega }} \right)}} }} - 2{\xi _0}} \right)\varepsilon _k^2{{\rm{\Omega }}_k} - Ch{\rm{\ }}{{\rm{\Omega }}_k}} \right]. \end{eqnarray*}$$

(22)

The effective mean stress ( |$\sigma _m^{\rm eff} = \ - \sigma _{kk}^{\rm eff}/3$|⁠) in (11b) was substituted here by the volumetric strain multiplied by the bulk modulus K ( |$\sigma _m^{\rm eff} = \ - K\ {I_1}$|⁠) leading to the power index N + 1.

The equation for the irreversible strain accumulation (11a) includes two terms. The first term, |$C_{ijnm}^{\psi \psi }\sigma _{nm}^{\rm eff}$|⁠, describes compaction/dilation with the rate proportional to the effective stress. The second is the damage-related coupling term multiplied to |${D_{ijkn}}$| (20), and describes the compaction or dilation associated with the formation and growth of damage. Representing the effective stress as a superposition of the volumetric (⁠|$\sigma _m^{\rm eff}$|⁠) and deviatoric (⁠|${\tau _{ij}}$|⁠) components allows us to describe the different mechanisms of the irreversible strain accumulation, or |${{\rm{\psi }}_{ij}}$| kinetics. The pressure driven compaction in the isotropic case becomes a well-known porosity reduction to its pressure-dependent equilibrium value, or Athy's (1930) law. Lyakhovsky et al. (2022) modified the scalar Athy relation and suggested that the 3-D equilibrium compaction depends on both pressure and deviatoric stress components:

$$\begin{eqnarray*} {\rm{\psi }}_{ij}^{\left( {\rm eq} \right)} = {B_0}\ \ \left[ {{\delta _{ij}} - \ \exp\left( { - \frac{{\sigma _m^{eff}}}{{{B_1}}}{\delta _{ij}} - \frac{{{\tau _{ij}}}}{{{B_2}}}} \right)} \right] \end{eqnarray*}$$

(23)

which further suggested the kinetics of the pressure-driven 3-D compaction has the form

$$\begin{eqnarray*} \frac{{{\rm d}{{\rm{\psi }}_{ij}}}}{{{\rm d}t}} = \ A\ \left( {{\rm{\psi }}_{ij}^{\left( {\rm eq} \right)} - {{\rm{\psi }}_{ij}}} \right)\ \sigma _m^{\rm eff}. \end{eqnarray*}$$

(24)

The eqs (23) and (24) consider not only the closure of voids or changes in the pore space (isotropic porosity reduction), but also changes in void shape under non-hydrostatic loading. Neglecting the term with deviatoric stress, or taking |${B_2} \to \infty $| in (23), reduces both the equilibrium compaction and kinetic equation to the traditional scalar form formulated in terms of material porosity.

Experimental studies suggest that permanent inelastic deformation is not only caused by pressure driven compaction, but also starts to accumulate at the onset of microcracking (as evidenced by the output of AE) and increases all the way up to the point of brittle failure (e.g. Lockner 1993b; Martin & Chandler 1994, 1998). This process is usually associated with the growth of microcracks and frictional sliding between grains, rather than closure of voids or space between grains. For similar reasons, Hamiel et al. (2004a) related the rate of irreversible strain accumulation with the rate of their scalar damage growth. Keeping in mind that the scalar damage variable is equivalent to the squared damage tensor, we extend their relation to:

$$\begin{eqnarray*} \frac{d}{{{\rm d}t}}\ \ \varepsilon _{ij}^{\left( {ir} \right)} = \left \{ \begin{array}{@{}*{1}{c}@{}} {{C_V}\frac{{d\left( {{{\rm{\Omega }}_{in}}{{\rm{\Omega }}_{jm}}} \right)}}{{{\rm d}t}}{\tau _{nm}},\ \ \ \ \frac{{d\left( {{{\rm{\Omega }}_{in}}{{\rm{\Omega }}_{jm}}} \right)}}{{{\rm d}t}} \gt 0}\\ {0,\ \ \ \ \ \ \ \ \ \ \ \frac{{d\left( {{{\rm{\Omega }}_{in}}{{\rm{\Omega }}_{jm}}} \right)}}{{{\rm d}t}} \le 0}. \end{array}\right. \end{eqnarray*}$$

(25)

The final kinetic eq. (11) for the |${\psi _{ij}}$| tensor incorporates all of the discussed mechanisms and to avoid a lengthy expression, is written for principal components:

$$\begin{eqnarray*} \frac{{{\rm d}{\psi _k}}}{{{\rm d}t}} = \ A\ \left( {{\rm{\psi }}_k^{\left( {\rm eq} \right)} - {{\rm{\psi }}_k}} \right)\ \sigma _m^{\rm eff} + {C_V}\frac{{d\left( {\Omega _k^2} \right)}}{{{\rm d}t}}{\tau _k} + {D_1}\exp\left[ { - {D_2}{{\rm{\Omega }}_k}} \right]\ {\left( { - {I_1}} \right)^N}\sqrt {{I_2}} \ \left[ {{\varepsilon _k}\sqrt {I_2^{\left( {\rm{\Omega }} \right)}} + \left( {\frac{{I_1^{\left( {\rm{\Omega }} \right)}}}{{\sqrt {I_2^{\left( {\rm{\Omega }} \right)}} }} - 2{\xi _0}} \right)\varepsilon _k^2{{\rm{\Omega }}_k} - Ch{\rm{\ }}{{\rm{\Omega }}_k}} \right]. \end{eqnarray*}$$

(26)

The first term of (26) represents the compaction prior to the onset of damage accumulation, which is a 3-D extension of the scalar Athy's compaction law. According to this term, the compaction approaches its stress-dependent equilibrium value with the rate proportional to the effective pressure. The second term is the 3-D equivalent of the damage-dependent irreversible strain accumulation with inverse of the effective viscosity or fluidity proportional to the rate of damage accumulation. This term describes extension or compaction depending on the sign of the deviatoric stress component. The last term represents the coupling between damage and porosity kinetics. Its sign, extension or compaction, is defined by the expression in the square brackets and depends on the loading and damage values.

The kinetic expressions for damage (22) and irreversible strain (26) provide the closed system of equations defining the 3-D evolution of the material properties.

2.6 Yield cap evolution

The new anisotropic poroelasic damage model shares the main features of the previously developed isotropic model of Lyakhovsky et al. (2015), including the ability to simulate the entire yield curve through a single formulation. Their model addressed several different deformation regimes including elastic deformation and pressure-driven compaction, brittle failure, and cataclastic flow. Loading of a rock sample to a level of stress beyond the initial yield surface caused an accumulation of damage and resulted in a porosity change (causing either compaction or dilation). In this case, the modelled yield cap grows. Consequently, if the sample is unloaded and then reloaded, the new yield cap is found to occur at a higher stress state, which is in agreement with the Kaiser effect. These general features of the model yield cap are shown in Fig. 1 using strain space, that is differential strain (⁠|${\varepsilon _1} - {\varepsilon _3}$|⁠) versus volumetric strain , instead more common stress space representation (differential stress versus pressure). Several authors formulated yield conditions in terms of strains and demonstrated significant advantages of this approach for materials with evolving yield conditions as a function of material properties and loading history (Naghdi & Trapp 1975; Yoder & Iwan 1981; Han & Chen 1986; Puzrin & Houlsby 2001). Although there are some similarities between stress-space and strain-space formulations, they are not equivalent when material weakening is considered (Casey & Naghdi 1983; Einav 2004, 2005). The yield surface (heavy black line in Fig. 1) represents the onset of the damage accumulation according to kinetic eq. (22) reduced to the model version for the isotropic material (⁠|${{\rm{\Omega }}_{ij}} = \ {\rm{\Omega }}\ {\delta _{ij}}$|⁠) and non-cohesive material (Ch = 0):

$$\begin{eqnarray*} \frac{{{\rm d}\alpha }}{{{\rm d}t}} = {C_d}\ \ \left[ {{D_1}{e^{ - {D_2}\alpha }}{{\left( { - {I_1}} \right)}^{N + 1}}\sqrt {{I_2}} + {I_2}\left( {\xi - {\xi _0}} \right)} \right] \end{eqnarray*}$$

(27)

Figure 1.

Three different deformational regimes and growing yield surface neglecting the cohesive force.

This form of the damage kinetic equation predicts the onset of damage accumulation (⁠|$d\alpha /{\rm d}t = {\rm{\ }}0$|⁠), or the yield condition expressed in terms of strain invariants:

$$\begin{eqnarray*} \sqrt {{I_2}} = \left[ {{D_1}{e^{ - {D_2}\alpha }}{{\left( { - {I_1}} \right)}^N} + {I_1}} \right]\ \ /\ {\xi _0}. \end{eqnarray*}$$

(28)

The entire yield curve (heavy black line in Fig. 1) is calculated for compactive volumetric strain (⁠|${I_1}$|⁠) between zero and a certain critical value |$I_1^*$|⁠, corresponding to the onset of damage under hydrostatic compaction. The differential strain (⁠|${\varepsilon _1} - {\varepsilon _3}$|⁠) were calculated using strain invariants and assuming triaxial loading conditions.The critical value (⁠|$I_1^*)$| is defined from the damage onset or yield condition (28) under hydrostatic loading corresponding to the strain invariant ratio |${I_1}/\ \sqrt {{I_2}} = \ - \sqrt 3 $|⁠:

$$\begin{eqnarray*} I_1^* = {\rm{\ }} - {\left( {\frac{{\sqrt 3 + {\xi _0}}}{{\sqrt 3 \cdot {D_1}{e^{ - {D_2}\alpha }}}}} \right)^{1/N}}. \end{eqnarray*}$$

(29)

The detailed discussion of the size and shape of the yield envelope and the sensitivity to the model parameters, including material cohesion, is presented in Appendix  B. The red line in Fig. 1 schematically represents the proportional load path with constant strain invariant ratio (⁠|$\xi \ = {I_1}/\sqrt {{I_2}} = \ {\rm Const.}$|⁠). At the initial stage of loading, when the stress level is beneath the yield cap (Regime I), quasi-elastic deformation is accompanied by material strengthening associated with crack closure and compaction (porosity decrease). When the load reaches the yield condition (red star), the damage accumulation starts. Distributed microcracking and grain crushing enable sliding along newly created internal surfaces, collapse of the pore space, and changes in grain packing arrangements leading to an overall porosity decrease. This deformational Regime II is most prominent under high confining pressures and is usually treated as cataclastic flow associated with compaction. Damage accumulation as well as porosity reduction leads to the decrease of the coupling coefficient D, which in the presented version is considered as damage-dependent only (eq. 8). The damage kinetics affects the yield condition (28) and the yield cap evolves as strain is accumulated. This feature allows the Kaiser effect to be reproduced, such that the onset of damage (and its associated AE) occurs at increasingly higher stress levels, if the sample is unloaded and reloaded along the same loading path.

Under elevated values of the differential strain, the loading path crosses the compaction-dilation transition (dashed line in Fig. 1) separating the deformational Regime II with compaction and Regime III with inelastic porosity increase (dilatancy). This transition coincides with the Coulomb–Mohr failure criterion and meets the condition |$\xi \ = {\xi _0}\ $|⁠. When the loading is pushed beyond the compaction–dilation transition, intensive damage accumulation, along with significant differential strains, lead to material dilation (porosity increase). Damage increase is then bounded by a certain critical value which eventually corresponds with macroscopic brittle failure, a dynamic stress drop and a rapid conversion of the differential elastic strain into plastic strain components.

The damage kinetic eq. (22) of the anisotropic model shares a similar structure to the scalar model (27). It also consists of two competing terms, but their values depend not only on the strain invariants, but also on the direction of loading relative to the principal values of the damage tensor. We rewrite these terms assuming that the principal directions of the damage tensor (⁠|${{\rm{\Omega }}_k}$|⁠) match the orientation of the principal loading axes with principal strain values |${\varepsilon _k}$|⁠:

$$\begin{eqnarray*} ter{m_1}\ \ = {D_1}\ \ \exp\left[ { - {D_2}{{\rm{\Omega }}_k}} \right]\ {\left( { - {I_1}} \right)^{N + 1}}\sqrt {{I_2}} \end{eqnarray*}$$

(30)

$$\begin{eqnarray*} ter{m_2} = \frac{1}{{\sqrt {{{\rm{\Omega }}_l}{{\rm{\Omega }}_l}} }}\ \ \left[ {{\varepsilon _k}\sqrt {I_2^{\left( {\rm{\Omega }} \right)}} + \left( {\frac{{I_1^{\left( {\rm{\Omega }} \right)}}}{{\sqrt {I_2^{\left( {\rm{\Omega }} \right)}} }} - 2{\xi _0}} \right)\varepsilon _k^2{{\rm{\Omega }}_k} - Ch{\rm{\ }}{{\rm{\Omega }}_k}} \right]. \end{eqnarray*}$$

(31)

We take out two constants, K and L, and rescale D₁ to preserve the same proportion between these terms. Fig. 2 demonstrates the evolution of these terms and their sums for the conventional loading, where |${{\rm{\Omega }}_1}$| is oriented in the direction of the axial load and |${\varepsilon _1}$| is the axial strain in the same direction. In this case, |${{\rm{\Omega }}_2},\ {{\rm{\Omega }}_3}$| are oriented in the direction of the transverse loads with two equal values of strain |${\varepsilon _2} = {\varepsilon _3}\ $|⁠. For the initially isotropic material ( |${{\rm{\Omega }}_1} = {{\rm{\Omega }}_2}\ \ = {{\rm{\Omega }}_3}\ $|⁠), the first term (30) is isotropic, since it depends on the damage and strain invariants. Its value strongly increases with the volumetric deformation (confinement) and weakly depends on the differential strain (Fig. 2a). Values of the second term (31) are however different depending on whether the axial or transverse directions are observed. This is the case even for the isotropic material (Figs 2 b1 and 2b2 since the value explicitly depends on the strain components. The axial values (Fig. 2) are always negative and increase (i.e. become more negative) with both volumetric and differential components. Thus, the sum of the first (30) and second (31) terms for the axial direction is negative for most of the strain values except for a small area in the right bottom corner of the map Fig. 2(c1). This implies that damage is not likely to occur in the axial direction. For most of the loading cases, except under high volumetric stress (strain) above critical value, |$I_1^*$| (29), the axial damage component remains unchanged or even decreases. Mechanically, this implies that microcracks which are oriented normal to the axial load direction become closed.

Figure 2.

Damage kinetics for initially isotropic damage (Ω₁ = Ω₂ = Ω₃) under conventional triaxial loading. (a) The 1st term is the same for all damage components. The 2nd term for axial component (b1) significantly differs from the values for the transversal components (b2). The damage kinetics in the axial direction (c1) differs from the transversal components (c2) which is similar in shape but not necessarily in values to the conventional scalar criterion.

The transverse values of the second term (31) are positive at relatively low volumetric strains and elevated differential strains (upper left-hand corner of Fig. 2b1). They become negative with volumetric strain increase. This change is indicated by a heavy red line separating the negative and positive values in Fig. 2(b2). The summation of the positive and negative values of the first (Fig. 2a) and second (Fig. 2b2) terms forms the yield cap (black heavy line in Fig. 2c2) which separates areas corresponding to damage accumulation (positive values) and healing or microcrack closure (negative values) for the transverse components |${{\rm{\Omega }}_2},\ {{\rm{\Omega }}_3}$|⁠.

The obtained shape of the yield cap for transverse damage components (red line in Fig. 3c2) is very similar to those predicted by the scalar model (see Figs B1-3 in Appendix  B). At low volumetric strains (with low level confinement) the yield cap is slightly shifted up depending on the cohesive force values (see Fig. B3 in Appendix  B for effect of cohesion). Loading above the yield cap leads to gradual accumulation of two transverse damage components, |${{\rm{\Omega }}_2},\ \ {{\rm{\Omega }}_3}$|⁠, that form at the same rate while the axial damage value, |${{\rm{\Omega }}_1}$|⁠, remains unchanged. This situation creates a stress-induced cylindrical transverse damage anisotropy observed in triaxial experiments (Browning et al. 2018). The damage can be schematically shown as two families of microcracks (insert in Fig. 3). One family (blue cracks in Fig. 3) is oriented parallel to the axial loading direction and opening in the transverse directions. Another family (black cracks in Fig. 3) normally oriented to the axial loading direction is closed. The axial damage value, |${{\rm{\Omega }}_1}$|⁠, remains unchanged until large hydrostatic volumetric strains above certain critical value, |$I_1^*$| close to 0.8 per cent, are applied. This value is similar to that predicted by the scalar model (29) and describes the onset of pore space collapse. The anisotropic model predicts that the axial damage increases under the loading conditions in the area on the right of the black dashed line in Fig. 3.

Figure 3.

Yield envelope for initially isotropic material.

Accumulation of the transverse damage components increases the size of the yield cap which is related to exponential damage-dependence in the first term (30). If the sample is then reloaded, along the same loading path, the damage commences at the point of the previous maximum stress value. This corresponds to the previously mentioned enlargement of the yield cap. This Kaiser effect is schematically shown in Fig. 4 as a ‘no rotation’ cyclic loading. During the first loading cycle (blue path and envelope), the onset of damage occurs at relatively low differential strain (strain values and model parameters will be specified in the next section comparing with experimental results). During the second cycle (purple path and envelope), the onset of damage occurs at a significantly higher value of strain (or stress), close to the maximum level of stress in the previous cycle. The quality of the Kaiser effect then depends on how close the yield cap keeps pace with the strain accumulation. This feature dramatically changes if the loading direction is rotated between consecutive loading cycles. Only a small change between the blue, green, and red yield envelopes is predicted for the case ‘rotation of the loading direction between cycles’. Similar to the previous case, the starting material has an essentially isotropic distribution of microcracks and flaws at the beginning of the first loading cycle (blue path and envelope in Fig. 4). However, after unloading at the end of the first cycle, the sample is no longer isotropic since new anisotropic microcrack damage has been formed, with the cracks growing parallel to the maximum loading direction (⁠|${{\rm{\Omega }}_2}$|⁠, |${{\rm{\Omega }}_3}$|⁠) and closure in the plane normal to the maximum loading direction (⁠|${{\rm{\Omega }}_1}$|⁠,). For the second loading cycle (red line), the loading orientations are rotated on the same sample that contains the previously formed anisotropic damage, such that the maximum loading direction is parallel to the previous |${{\rm{\Omega }}_2}$|⁠,. The low Ω₁ and one of the high Ω₃ damage components are now associated with the transverse direction, while another high component Ω₂ is axial. The onset of accumulation of the new Ω₁ damage component occurs at almost exactly the same level of stress as in the first cycle. There is no observation of any Kaiser effect, and this suggests that a completely new family of microcracks are generated normal to those generated on during the first cycle. An exactly similar scenario is observed during the third loading cycle when the sample is reloaded with the loading orientation again rotated such that the |${{\rm{\Omega }}_3}$| component becomes parallel to the new maximum loading direction (the green path). The onset of damage again occurs at essentially the same level of stress and strain as in cycle 1 and cycle 2. Again, no Kaiser effect is observed. Analysis of Fig. 4 explains why there is apparently no Kaiser effect. The blue, red and green lines related to the first, second and third sequential cycles exhibit very similar yield caps which suggests that the new damage on each cycle is independent of the damage formed during earlier cycles. We would therefore not expect to encounter a Kaiser effect under these conditions. It is only when the sample is loaded to a higher level of stress in the same orientation, as demonstrated by the purple line, that would expect to observe a manifestation of the Kaiser effect.

Figure 4.

Evolving yield envelope for loading cycles in which the maximum load is applied in the same direction cyclically and in which the maximum loading direction is sequentially rotated with respect to the sample.

To explain this model prediction further, Fig. 5 shows the distribution of damage kinetics terms for the two transverse low Ω₁ and high Ω₃ damage components during the second load cycle for the sample with anisotropic damage Ω₁ < Ω₂ = Ω₃. In spite of the elevated damage in the axial direction, similar to the isotropic sample, the axial damage component, Ω₂, remains unchanged or even decreased or healed (maps are not shown). Values of the first (Figs 5a1 and a2) and the second (Figs 5b1 and b2) terms (30, 31) are similar. However, the first term for the Ω₁ component (Fig. 5a2) is slightly larger than for the Ω₃ (Fig. 5a1). The values of the second term show the opposite tendency (Figs 5b1 and b2). Finally, the sum of these terms controlling the damage kinetics (Fig. 5c1 and c2), gives significantly different yield cap (red heavy lines in Figs 5c1 and c2). Fig. 6 summarizes the shape and size of the yield envelope for the anisotropic material (red lines) and compares them with yield caps for the initially isotropic case (dotted lines). The yield cap of the low transverse damage component (red line for Ω₁) is essentially coincident with the yield cap for the isotropic material. However, the yield cap for the high transverse damage component Ω₃ is significantly enlarged. If the loading is not large enough and it is only slightly above the red line for Ω₁ component, the high transverse component Ω₃ remains unchanged and only minor total damage is accumulated in the sample. The yield for the axial Ω₂ component is shifted to significantly higher values of volumetric strain, value, |$I_1^*$| above 1.0 per cent, but retaining the same shape as the yield for the axial component of the isotropic material.

Figure 5.

Kinetics of the transverse damage components (Ω₁, Ω₃) for initially anisotropic initial damage (Ω₁ < Ω₂ = Ω₃). Upper raw represents the values for the larger (Ω₃) component and the lower raw for the smaller (Ω₁) component. The pattern for both components, 1st (a1, a2) and 2nd (b1, b2) terms are similar, but their different values lead to significantly different size of the yield envelope (c1, c2).

Figure 6.

Yield envelope for initially anisotropic material (Ω₁ < Ω₂ = Ω₃). The envelope for the smallest transverse (Ω₁) component (red line) is almost the same as for the initially isotropic material (dotted blue line). The envelope for larger transversal component (Ω₃) is significantly larger than for the component Ω₁. The axial damage component (Ω₂) is accumulated only under high-confining conditions.

3 MODEL VERIFICATION

3.1 Materials and experimental settings

In order to verify the ability of our new model to reproduce the features of the experimentally observed 3-D Kaiser effect, we use results from two sets of tests reported by Browning et al. (2018), namely; the sequential rotational conventional triaxial (SRCT) loading test, and the cyclic SRCT (CSRCT) loading test. The experiments were conducted on dry samples of the relatively homogeneous Darley Dale sandstone, which is a feldspathic sandstone with a moderate porosity of ∼13 per cent and a grain size range from 0.08 to 0.8 mm (Wu et al. 2000; Heap et al. 2009). The deformation apparatus, based at the laboratories of Koninklijke Shell Exploratie en Produktie Laboratoriu (KSEPL), Rijswijk, Netherlands, consists of a three-axis stressing frame constructed of flanged steel beams, one of which was removable to allow the insertion of the cubic rock samples (edge length of 50 mm). Loading was performed with three pairs of servo-controlled hydraulic rams with a loading capacity of 300 kN. Hemispherical seatings were used along orthogonal axes perpendicular to the faces of the cubic samples and loading platens, with an edge length of 47.5 mm were interposed between the rams and the sample faces to provide the contact surfaces.

In order to keep the sample centered within the apparatus and to ensure good acoustic contact between the sample and the loading platens, a small pre-load of 4 MPa was applied along each of the three axes prior to testing. The load in each of the three directions was measured using electronic load cells with an accuracy of ±0.2 per cent, and the displacement in each direction was measured using linear variable differential transformers (LVDTs) mounted between the loading platens. AE was monitored using a piezo-electric transducer located in a recess within one of the platens.

3.2 SRCT loading

The SRCT loading test consisted of three loading cycles (Fig. 7a). Prior to the first loading cycle, the sample was pre-loaded hydrostatically up to 4 MPa. Then, in the first loading cycle, the differential stress (blue lines in Fig. 7a) was increased in the 1-direction to a maximum value of 80 MPa at a rate of 0.018 MPa s^–1 and then unloaded at the same rate, while both transverse stresses were held constant at 4 MPa. To relate the experimental results to our model formulation, the axial loading direction during the first cycle is associated with the first damage component Ω₁, while the transverse damage components are Ω₂ and Ω₃. In the second cycle, the differential stress was rotated to the 2-direction (red lines in Fig. 7a), and the new damage then corresponds to damage component Ω₂, while the transverse components become Ω₁ and Ω₃. The loading protocol was the same as in the first cycle such that the differential load was increased at a rate of 0.018 MPa s^–1 to a maximum value of 80 MPa and then unloaded at the same rate. Again, the two transverse stresses were held constant at 4 MPa. Finally, the same loading protocol was applied in the third cycle following a further rotation of the differential to the 3-direction (green lines in Fig. 7a), such that Ω₃ now corresponds to the axial component and Ω₁ and Ω₂ correspond to the transverse components.

Figure 7.

Comparison between experimental results and modelling. (a) Sequential rotational conventional triaxial loading; markers show the onset of the observed AE (black line in b). Coloured line in (b) shows the simulated damage accumulated during the loading. (c) Rate of damage accumulation

Fig. 7(b) (black line) demonstrates that the onset of AE occurs at approximately the same level of stress in each cycle; between 35 and 45 MPa. This suggests that no manifestation of the Kaiser effect is observed in this test. However, while the onset of AE occurs at approximately the same level of stress in each cycle, it is significant to note that the amount of AE (plotted as cumulative AE hits in Fig. 7b), decreases with each new sequential loading cycle. The model-simulated accumulated damage curves recreate this tendency for the three loading cycles (colored lines for each cycle in Fig. 7b). The simulation commences with an initially small amount of isotropic damage, |${{\rm{\Omega }}^2} = {\rm{\Omega }}_1^2\ \ + {\rm{\Omega }}_2^2 + {\rm{\Omega }}_3^2 \approx 1.5\ {\rm per\ cent}$|⁠. During the first loading cycle, the two transverse damage components, Ω₂ and Ω₃, grow at the same rate (as shown in Fig. 7c) leading to a damage increase, |${{\rm{\Omega }}^2} \ge 4\ {\rm per\ cent}$|⁠. However, during the second and third loading cycles, the two transverse damage components (Ω₁ and Ω₃ in cycle 2, Ω₁ and Ω₂ in cycle 3) grow at different rates. This is entirely as expected, because the rock is no longer isotropic following the anisotropic damage formed during the first cycle. The overall amount of damage increase during cycles 2 and 3 together add only slightly more than 1 per cent to the total amount of damage. Finally, the simulated stress–strain curves from the model are shown in Fig. 8 (coloured lines) and provide a very good fit to the experimental curves (black lines). However, it should be noted that the fitting procedure is done by eye and the search for model parameters is therefore non-unique. There are then several uncertainties and trade-offs between the model parameters. For example, the stresses and strains prior to the onset of damage accumulation (AE output) for the low damage starting material are expected to be linearly elastic, and hence permit the calculation of elastic moduli. However, for this specific test, on a relatively porous sandstone, a significant amount of the total strain is irreversible and associated with material compaction (eqs 23, 24). Panteleev et al. (2021), who processed results for similar experiments on samples of Darley Dale sandstone, also noted this feature. We therefore use Lamé moduli |${\lambda _0} = \ 5\ {\rm GPa},\ \ \ {\mu _0} = \ 7\ {\rm GPa}$| for the starting material. These values are in agreement with the measured seismic wave velocities of V_p ∼ 3.4 km s^–1 and V_s ∼ 2.1 km s^–1 reported by Browning et al. (2017) for this material. Since all tests were performed under the same hydrostatic confining pressure of 4 MPa, it is impossible to constrain the shape of the yield curve. With this uncertainty we use the power index, N = 1, in (20); |${\xi _0} = \ - 0.1$| and Ch = 2 10⁻⁷ in (22) which corresponds to about 15 MPa of cohesive force (Appendix  B). Scaling the range of the damage value, |${{\rm{\Omega }}^2}$|⁠, from zero to 100 per cent for total failure (e.g. Lyakhovsky et al. 1997) gives |$\gamma \ = \ 10\ {\rm GPa}$|⁠. Simultaneous fitting of the stresses and strains in the three cycles and the noted similarity between the accumulated damage and the cumulative AE output led to the following parameters: compaction equilibrium, B₀ = 0.5 per cent, B₁ = 20 MPa, B₂ = 15 MPa in (23), C_v = 3 × 10⁻² (MPa s)⁻¹ in (25), rate coefficient in (24) A = 10⁻⁵ (MPa s)⁻¹, together with L = 7 s⁻¹, K × D₁ = 7 × 10³ s⁻¹ and D₂ = 30 in (22). The relatively high D₂ value corresponds to the high sensitivity of the yield curve to the change in damage required to reproduce the Kaiser effect (e.g. Gajst et al. 2020).

Figure 8.

Comparison between observed (coloured lines) and simulated (grey lines) stress–strains for three loading cycles (blue, red, green).

3.3 CSRCT loading

CSRCT loading was performed using the same starting material and with the same initial pre-loading of 4 MPa along each of the three sample axes. In the first phase of this test, the differential stress was initially raised to 75 MPa along the 1-direction (i.e. the direction of the Ω₁ damage component) using the same loading rate as for the SRCT test, and the sample was then unloaded instantaneously (Fig. 9a). The sample was subsequently reloaded in the same orientation but to a higher differential stress of 80 MPa (blue lines in Fig. 9a), before again being unloaded instantaneously. Then, during the second phase, the differential stress was rotated and applied in the 2-direction (i.e. the direction of the Ω₂ damage component) to a level of 65 MPa in the first cycle and 80 MPa in the second cycle (red lines in Fig. 9a). Finally, in the third phase, the differential stress was rotated again and applied in the 3-direction (i.e. the direction of the Ω₃ damage component). During this phase, the sample was loaded to 55 MPa in the first cycle and 80 MPa in the second cycle (green lines in Fig. 9a). Once again, we see that the model simulated damage curves obtained using the same material properties as for the SRCT test, are very similar in form to the cumulative AE output recorded during the experiment (Fig. 9b). Similarly to the SRCT test, the two transverse damage components, Ω₂ and Ω₃, grow at the same rate (Fig. 9c) during the loading in the 1-direction, as expected for the isotropic sample. However, during the loading in 2- and 3-directions the rock is no longer isotropic following the anisotropic damage formed during the loading in 1-direction. After stress was rotated, the two transverse damage components Ω₁ and Ω₃, and then, after another rotation, components Ω₁ and Ω₂ grow at different rates. The modeled onsets of damage for each loading cycle, are in very good agreement with the measured onsets of AE output for each cycle during the experiment. The colored markers (stars) on the loading curves of Fig. 9(a) indicate the value of stress at the onset of AE, and hence the onset of new damage. We note that both the experimental data and the model simulations exhibit a distinct Kaiser effect when the samples were reloaded in the same direction, but no such effect when the differential stress was rotated. Again, this suggests that the onset of damage is, at best, only weakly affected by damage accumulated during earlier phases when the differential stress was applied in different orientations.

Figure 9.

Comparison between experimental results and modelling. (a) Cyclic sequential rotational conventional triaxial loading; the markers show the onset of the observed AE (black line in b). Coloured line in (b) shows the simulated damage accumulated during the loading. (c) Rate of damage accumulation.

4 DISCUSSION AND CONCLUDING REMARKS

Here we build on the previous non-linear anisotropic damage rheology model (Panteleev et al. 2021) by presenting a newly developed poroelastic rheological model which accounts for both coupled anisotropic damage and porosity evolution. The new model shares the main features of our previously developed anisotropic damage and scalar poroelastic damage models, including the ability to simulate the entire yield curve through a single formulation. In the new model, the yield condition is defined in terms of invariants of the strain tensor, and so the new formulation operates with directional yield conditions (different values for each principal direction) depending on the damage tensor and triaxial loading conditions. This allows us to discern evolving yield conditions for each principal stress direction and fit the measured amounts of accumulated damage from previous loading cycles. Coupling between anisotropic damage and anisotropic compaction along with the damage-dependent yield condition produces a reasonable fit to the experimentally obtained stress–strain curves. Furthermore, the simulated time-dependent cumulative damage is well correlated with experimentally observed AEs during cyclic loading in different directions. As such, we are able to recreate many of the features of the experimentally observed directional 3-D Kaiser ‘damage memory’ effect.

The main finding from this formulation is that each independent direction will posses its own yield envelope. The state of each envelope then depends on the direction, magnitude and history of loading. Yield in the transverse components are similar in shape but not necessarily in value to the conventional scalar criterion and the yield criterion in the axial direction has an entirely different shape. As such, damage accumulation in the axial direction is possible only under high volumetric stresses.

These results are important in nature since rocks in complex tectonic environments, such as volcanoes or active fault zones, experience stresses that evolve both spatially and temporally and experience not only cyclic loading and unloading, but also rotation and/or reorientation of stresses. The resulting crack distributions will then form sequentially and may be highly anisotropic. Thus, the tectonic history of the crust as recorded in deformed rocks may include evidence for complex stress paths, encompassing different magnitudes and orientations. Geodetic and seismic data from periods of inflation and deflation at Krafla volcano in Iceland demonstrates that the rate of seismicity increases only after the amount of inflation in a previous cycle has been reached or exceeded (Heimisson et al. 2015). These, and similar observations, point to a potential crustal scale Kaiser effect but the directionality of fracture populations formed has received less attention. The data from Heimisson et al. (2015) were interpreted assuming a conventional Kaiser effect such that the orientation of loading and unloading of all of the episodes were assumed to stress the crustal rocks in the same direction. It is noted that, prior to many cyclic inflation episodes at active volcanoes, the level of inflation has often been at a higher level the than the first cycle of a new episode. In such circumstances, the first cycle should not produce seismicity, as the previous level of inflation has not been reached or exceeded. This has traditionally been explained as related to crack healing processes between cycles (Kim et al. 2014) but it can also be explained by our new 3-D model. As the orientation and magnitude of crustal stresses can vary enormously between different rock layers (Gudmundsson 2011) the new inflation episode may have loaded the rocks under a slightly different axis and hence triggered seismicity at lower stresses, similar to as demonstrated in Fig. 9(a).

Recent experimental results have shown that damage, under true triaxial loading, is a distinctly directional phenomenon, and these results have also revealed a 3-D directionally dependent Kaiser ‘damage memory’ effect. The developments of this study provide an internally consistent framework for simulating evolving crustal rock damage under repeated cycles of stress in complex tectonic environments where stresses may evolve both spatially and temporally.

ACKNOWLEDGEMENTS

The paper benefited from useful comments by two referees, Manolis Veveakis and Klaus Regenauer-Lieb, and the editor, Alexis Maineult. The contributions by Lyakhovsky and Shalev was supported by grant from the Israel Science Foundation, ISF 363/20. The contributions by Browning, Meredith, Healy and Mitchell were supported by UKRI NERC awards NE/N003063/1, NE/N002938/1, NE/T007826/1 and NE/T00780X/1. The contributions by Browning was also supported by FONDECYTgrant number 11190143. The contribution by Panteleev was supported by Russian Science Foundation (project N 19-77-30008).

DATA AVAILABILITY

Experimental data may be obtained from J.B. (e-mail: [email protected]).

APPENDIX A

Thermodynamic relations

The presented formulation is based on balance equations of the irreversible thermodynamics of the continuum media with internal variables (Malvern 1969; Coussy 1995). We consider the system with density of the free energy, F, to be a function of:

$$\begin{eqnarray*} F\ = \ F\left( {T,{{\rm{\varepsilon }}_{ij}},{{\rm{\Omega }}_{ij}},{{\rm{\psi }}_{ij}},{\rm{\zeta }}} \right), \end{eqnarray*}$$

(1)

where T is the temperature, |${{\rm{\varepsilon }}_{ij}}$| is the elastic strain tensor, |${{\rm{\psi }}_{ij}}$| is the compaction tensor, |${{\rm{\Omega }}_{ij}}$| is the damage tensor and |${\rm{\zeta }}$| is the change in volume fluid content defined by Biot (1941), see also Detournay & Cheng (1993). Since each variable can vary independently of the other variables, Gibbs relation can be written as (Gibbs 1961):

$$\begin{eqnarray*} dF\ = \ - SdT + \frac{{\partial F}}{{\partial {\varepsilon _{ij}}}}{\rm d}{\varepsilon _{ij}} + \frac{{\partial F}}{{\partial {\Omega _{ij}}}}{\rm d}{\Omega _{ij}} + \frac{{\partial F}}{{\partial {\psi _{ij}}}}{\rm d}{\psi _{ij}} + \frac{{\partial F}}{{\partial \zeta }}{\rm d}\zeta , \end{eqnarray*}$$

(2)

where |$S\ = \ - \frac{{\partial F}}{{\partial T}}$| is entropy density (Einstein's summation convention is assumed).

The balance equation for the density of the internal energy, U, includes three source terms associated with shear heating (stress, |${{\rm{\sigma }}_{ij}}$|⁠, times strain rate or time derivative of the total strain rate tensor, |${{\rm{e}}_{ij\ }}$|⁠), energy dissipation due to heat flux, |${Q_i}$|⁠, and the advective flux due to fluid flow (eq. 71 from Coussy et al. 1998):

$$\begin{eqnarray*} \frac{{dU}}{{{\rm d}t}} = \frac{d}{{{\rm d}t}}\ \ \ \ \left( {F + TS} \right) = {{\rm{\sigma }}_{ij}}\ \ {{\rm{e}}_{ij\ }} - {\nabla _i}{Q_i} - {\nabla _i}\left[ {\ \left( {{h_f} + T{s_f}} \right){q_i}} \right], \end{eqnarray*}$$

(3)

where |${q_i}$| is fluid flux, |${h_f}$| is the enthalpy or Gibbs potential of the fluid, which is a function of the fluid pressure, |${p_f}$|⁠, and temperature; |${s_f}$| is the entropy of the fluid. Similarly, the entropy balance equation includes positive entropy production, |$\Gamma $|⁠, divergence of the heat flux, and similar advective term due to fluid flow (eq. 72 from Coussy et al. 1998):

$$\begin{eqnarray*} \frac{{dS}}{{{\rm d}t}} = \frac{\Gamma }{T}\ \ - {\nabla _i}\left( {\frac{{{Q_i}}}{T}} \right) - {\nabla _i}\left( {{q_i}\ {s_f}} \right). \end{eqnarray*}$$

(4)

The stress tensor and the fluid pressure are defined as (Malvern 1969; Coussy 1995):

$$\begin{eqnarray*} {{\rm{\sigma }}_{ij}} = \frac{{\partial F}}{{\partial {\varepsilon _{ij}}}}{\rm{\ }} \end{eqnarray*}$$

(5)

$$\begin{eqnarray*} {p_f} = \frac{{\partial F}}{{\partial \zeta }}{\rm{\ }}. \end{eqnarray*}$$

(6)

Fluid mass conservation is:

$$\begin{eqnarray*} \frac{{{\rm d}\zeta }}{{{\rm d}t}} = \ - \ {\nabla _i}\left( {{q_i}} \right). \end{eqnarray*}$$

(7)

Combining (2, 3, 4) and using (5, 6, 7) we get a final expression for the local entropy production in a form:

$$\begin{eqnarray*} \Gamma \ = \ - \frac{{{Q_i}}}{T}{\nabla _i}T - {q_i}{\nabla _i}{p_f} + {{\rm{\sigma }}_{ij}}\left( {{{\rm{e}}_{ij\ }} - \frac{{{\rm d}{\varepsilon _{ij}}}}{{{\rm d}t}}} \right) - \frac{{\partial F}}{{\partial {\Omega _{ij}}}}\frac{{{\rm d}{\Omega _{ij}}}}{{{\rm d}t}} - \frac{{\partial F}}{{\partial {\psi _{ij}}}}\frac{{{\rm d}{\psi _{ij}}}}{{{\rm d}t}}. \end{eqnarray*}$$

(8)

First and second terms represent dissipation associated with heat transport and fluid flow. The Fourier and Darcy laws establish linear relations between temperature and fluid pressure gradients with corresponding fluxes:

$$\begin{eqnarray*} {Q_i}\ \ = - K_{ij}^T{\nabla _j}\ \ T \end{eqnarray*}$$

(9a)

and

$$\begin{eqnarray*} {q_i}\ \ = - K_{ij}^F{\nabla _j}\ \ {p_f}, \end{eqnarray*}$$

(9b)

where |$K_{ij}^T$| and |$K_{ij}^F$| are positively defined thermal conductivity and permeability tensors.

The elastic strain tensor is the difference between total, |${g_{ij}}$|⁠, and irreversible compaction, |${{\rm{\psi }}_{ij}}$|⁠, tensors:

$$\begin{eqnarray*} {{\rm{\varepsilon }}_{ij}} = {g_{ij}}\ \ - {{\rm{\psi }}_{ij}}. \end{eqnarray*}$$

(10)

Taking the time derivative of (10) and substituting the total strain rate tensor, |${{\rm{e}}_{ij}} = {\rm d}{g_{ij}}/{\rm d}t\ ,$| into (8), the part of the total dissipation, |${\Gamma _{DC}}$|⁠, associated with evolving damage and compaction (three last term of eq. 8) is expressed as a sum of two terms proportional to the damage and compaction rates:

$$\begin{eqnarray*} {\Gamma _{DC}} = \ - \frac{{\partial F}}{{\partial {\Omega _{ij}}}}\frac{{{\rm d}{\Omega _{ij}}}}{{{\rm d}t}} + \left( {{{\rm{\sigma }}_{ij}} - \frac{{\partial F}}{{\partial {\psi _{ij}}}}} \right)\frac{{{\rm d}{\psi _{ij}}}}{{{\rm d}t}} > 0. \end{eqnarray*}$$

(11)

For small deviations from equilibrium, the entropy production or dissipation potential may be approximated as a quadratic function of the rate of the internal variables. In this case constitutive equations give the thermodynamic fluxes as a linear function of the thermodynamic forces (Malvern 1969; DeGroot & Mazur 2013) also known as the Onsager's (1931) relations:

$$\begin{eqnarray*} \frac{{{\rm d}{\psi _{ij}}}}{{{\rm d}t}} = C_{ijnm}^{\psi \psi }\ \ \left( {{{\rm{\sigma }}_{ij}} - \frac{{\partial F}}{{\partial {\psi _{ij}}}}} \right) - C_{ijnm}^{\psi \Omega }\frac{{\partial F}}{{\partial {\Omega _{ij}}}} \end{eqnarray*}$$

(12a)

$$\begin{eqnarray*} \frac{{{\rm d}{\Omega _{nm}}}}{{{\rm d}t}} = C_{ijnm}^{\Omega \psi }\ \ \left( {{{\rm{\sigma }}_{ij}} - \frac{{\partial F}}{{\partial {\psi _{ij}}}}} \right) - C_{ijnm}^{\Omega \Omega }\frac{{\partial F}}{{\partial {\Omega _{ij}}}}. \end{eqnarray*}$$

(12b)

These phenomenological kinetic equations guarantee the non-negative entropy production if the cells of the matrix of the kinetic coefficients

$$\begin{eqnarray*} {C_{ijnm}} = \left| {\begin{array}{@{}*{2}{c}@{}} {C_{ijnm}^{\psi \psi }}&{C_{ijnm}^{\psi \Omega }}\\ {C_{ijnm}^{\Omega \psi }}&{C_{ijnm}^{\Omega \Omega }} \end{array}} \right|{\rm{\ }} \end{eqnarray*}$$

(13)

meet several conditions (Malvern 1969; DeGroot & Mazur 2013). Matrixes of the diagonal cells (⁠|$C_{ijnm}^{\psi \psi }$|⁠, |$C_{ijnm}^{\Omega \Omega }$|⁠) must be positively defined. Off-diagonal terms are usually taken to be either symmetric or antisymmetric. Following poroelastic damage model of Hamiel et al. (2004b) and Lyakhovsky et al. (2015) we adopt antisymmetric structure ( |$C_{ijnm}^{\Omega \psi } = \ - C_{ijnm}^{\psi \Omega }$|⁠) of the kinetic matrix (13). These conditions assure positive dissipation (11).

Further model formulation, demonstrating its main features as well as its calibration, verification and application to specific problems require definition of the energy function (1) and kinetic coefficients (13). In the Appendix  B we briefly discuss the isotropic, scalar formulation, which instead of the tensorial damage-compaction variables operates with the scalar damage, |$\alpha $|⁠, equal to the squared value, |$\alpha \ = \ {{\rm{\Omega }}^2}$|⁠, derived by Panteleev et al. (2021), and porosity, |$\varphi $|⁠. Scalar |${\rm{\Omega }}$| and |$\varphi $| variables are connected with the tensor variables as:

$$\begin{eqnarray*} {{\rm{\Omega }}_{ij}} = \ {\rm{\Omega }}\ {\delta _{ij}} \end{eqnarray*}$$

(14a)

$$\begin{eqnarray*} {{\rm{\psi }}_{ij}} = \ \varphi \ {\delta _{ij}}. \end{eqnarray*}$$

(14b)

The complete anisotropic formulation is presented in the Appendix  C.

APPENDIX B

Scalar poroelastic damage model

Following Biot's theory of poroelasticity (Biot 1941, 1956) the free energy of a poroelastic medium, F, is a sum of the elastic energy and the poroelastic coupling term of the saturated medium with the Biot modulus, M, and the Biot coefficient for porous media, |$\beta $|⁠:

$$\begin{eqnarray*} F\ = \frac{{\lambda \left( {\alpha ,\varphi } \right)}}{2}\ \ I_1^2 + \mu \left( {\alpha ,\varphi } \right){I_2} - \gamma \left( {\alpha ,\varphi } \right){I_1}\sqrt {{I_2}} + \frac{1}{2}M{\left[ {\beta {I_1} - \zeta + \varphi } \right]^2} + Ch\ {\gamma _m} \alpha . \end{eqnarray*}$$

(15)

In the literature discussing Biot poroelasticity, this coefficient is often noted as |${\rm{\alpha }}$|⁠. To avoid a duplicate notation we use |$\alpha $| for the scalar damage and change the notation for the Biot coefficient to |$\beta $|⁠.

Since the target model is isotropic, the energy function may depend only on invariants of the elastic strain tensor, |${{\rm{\varepsilon }}_{ij}}$| ( |${I_1} = {{\rm{\varepsilon }}_{ii}}\ ,\ \ \ \ {I_2} = {{\rm{\varepsilon }}_{ij}}\ \ {{\rm{\varepsilon }}_{ij}}$|⁠). Following (Hamiel et al. 2004b, Lyakhovsky et al. 2015) the elastic energy for nonlinear damaged media includes two Hookean terms with the Lame drained moduli λ, μ  and an additional non-linear term with strain coupling modulus γ. Value of this additional modulus varies from zero for damage free material with |$\alpha \ = \ 0$| to |$\gamma \ = {\gamma _m}\ \ $| at material failure with |$\alpha \ = \ 1$| (Lyakhovsky et al. 1997). Following Gajst (2020), we introduce the damage-dependent term with the non-dimensional coefficient Ch. Below we will show that this term allows accounting for the cohesive force which is important for loading conditions under low confining pressures.

Differentiation of the poroelastic energy (15) according to the constitutive relations (5, 6), the stress tensor, |${{\rm{\sigma }}_{{\rm{ij}}}}$|⁠, and fluid pressure, |${{\rm{p}}_{\rm{f}}}$|⁠, are:

$$\begin{eqnarray*} {{\rm{\sigma }}_{ij}} = \frac{{\partial F}}{{\partial {\varepsilon _{ij}}}}\ \ = \left( {\lambda - \frac{\gamma }{\xi }} \right)\ \ \ {I_1}{{\rm{\delta }}_{ij}} + \left( {2\mu - \gamma \xi } \right){{\rm{\varepsilon }}_{ij}} + \beta M\left( {\beta {I_1} - \zeta + \varphi } \right) \end{eqnarray*}$$

(16)

$$\begin{eqnarray*} {p_f} = \frac{{\partial F}}{{\partial \zeta }}\ \ = \ M\left( { - \beta {I_1} + \zeta - \varphi } \right), \end{eqnarray*}$$

(17)

where |$\xi \ = {I_1}\ /\sqrt {{I_2}} $| is the strain invariant ratio changing from |$\xi = - \sqrt 3 $| for isotropic compaction to |$\xi = \sqrt 3 $| for isotropic dilation. Using the assumption that |$\lambda \ = \ Const.$| and both |$\mu ,\ \gamma $| linearly depend on the damage (Agnon & Lyakhovsky 1995), the derivatives of the energy function (15) are the effective stress:

$$\begin{eqnarray*} {\sigma _{ij}} - \frac{{\partial F}}{{\partial \varphi }}\ \ \ {\delta _{ij}} = \sigma _{ij}^{\rm eff}{\rm{\ }} \end{eqnarray*}$$

(18)

and damage induced energy change:

$$\begin{eqnarray*} \frac{{\partial F}}{{\partial \alpha }} = {\gamma _m}\ \ \left( { - \ {I_2}\ \left( {\xi - {\xi _0}} \right) + Ch} \right), \end{eqnarray*}$$

(19)

where |${\xi _0}$| is the critical value of the strain invariant ratio corresponding to the conditions of the Coulomb failure criteria (see Lyakhovsky et al. 1997, for details). Substituting the derivatives (18, 19) into kinetic eq. (12) and accounting to the scalar nature of the damage, |$\alpha $|⁠, and porosity, |$\varphi $|⁠, variables (14), and using the mean stress|$( {\ \ {\sigma _m} = \ - {\sigma _{kk}}/3} )$|⁠, the coupled kinetic equations are reduced to (see Hamiel et al. 2004b, Lyakhovsky et al. 2015, and references therein):

$$\begin{eqnarray*} \frac{{{\rm d}\varphi }}{{{\rm d}t}} = \ - A\ {\sigma _m} - D\ \left[ {{I_2}\left( {\xi - {\xi _0}} \right) - Ch} \right] \end{eqnarray*}$$

(20)

$$\begin{eqnarray*} \frac{{d\alpha }}{{{\rm d}t}} = \ D\ {\sigma _m} + {C_d}\left[ {{I_2}\left( {\xi - {\xi _0}} \right) - Ch} \right]. \end{eqnarray*}$$

(21)

Hamiel et al. (2004b) suggested that the coupling coefficients D is a power-law expression of the effective pressure, |$D\ \sim \ \sigma _m^N$|⁠. They demonstrated that the transition from positive to negative values of the slope of the yield curve (yield cap) is a general feature of the model. Here we modify the previous formulation using strain invariants |${I_1},\ \ \sqrt {{I_2}} $| and |$\ \ \ {\sigma _m} = \ - K\ {I_1}$| (K is a bulk modulus); and also following Gajst et al. (2020) D-value exponentially decreases with damage:

$$\begin{eqnarray*} D\ = {C_d}\ \ {D_1}\ \ {e^{ - {D_2}\alpha }}\ {\left( { - {I_1}} \right)^N}\sqrt {{I_2}} /K. \end{eqnarray*}$$

(22)

N > 0, D₁ and D₂ values define the shape of the yield envelope. The form (22) is defined only for compaction (⁠|${I_1} \le 0$|⁠); D₁ is either constant or function of the porosity constrained using strain-defined yielding envelope. With this assumption, the kinetic for damage accumulation (21) is:

$$\begin{eqnarray*} \frac{{d\alpha }}{{{\rm d}t}} = {C_d}\ \ \left[ {{D_1}{e^{ - {D_2}\alpha }}{{\left( { - {I_1}} \right)}^{N + 1}}\sqrt {{I_2}} + {I_2}\left( {\xi - {\xi _0}} \right) - Ch} \right]. \end{eqnarray*}$$

(23)

This form of the damage kinetic equation predicts the onset of damage accumulation ( |$\frac{{d\alpha }}{{{\rm d}t}} = {\rm{\ }}0$|⁠) for the intact material (⁠|$\alpha \ = \ 0$|⁠), or the yield condition:

$$\begin{eqnarray*} {D_1}{\left( { - {I_1}} \right)^{N + 1}}\sqrt {{I_2}} + {I_2}\left( {\xi - {\xi _0}} \right) - Ch\ = {\rm{\ }}0. \end{eqnarray*}$$

(24)

Under shear load needed to overcome the cohesive force ( |${\tau _h} = \ 2{\mu _0}{\varepsilon _h}$|⁠) with zero volumetric strain ( |${I_1} = \ 0,\ \ \xi \ = \ 0\ {\rm{and}}\ \ \ {I_2} = \ 2\varepsilon _h^2$|⁠), the coefficient Ch is equal to:

$$\begin{eqnarray*} Ch\ = {\rm{\ }} - {\xi _0}\frac{{\tau _h^2}}{{2\mu _0^2}}. \end{eqnarray*}$$

(25a)

Under large hydrostatic volumetric strains, neglecting the cohesive force, the onset of damage occurs at:

$$\begin{eqnarray*} I_1^* = {\rm{\ }} - {\left( {\frac{{\sqrt 3 + {\xi _0}}}{{\sqrt 3 \cdot {D_1}}}} \right)^{1/N}}. \end{eqnarray*}$$

(25b)

The eqs (25a) and (25b) define two endmember yield values. The entire yield shape may be calculated for any given volumetric strain between zero and |$I_1^*$| by solving (24) as a quadratic equation for |$\sqrt {{I_2}} $|⁠:

$$\begin{eqnarray*} \sqrt {{I_2}} = \frac{{{D_1}{{\left( { - {I_1}} \right)}^N} + {I_1} + \sqrt {{{\left( {{D_1}{{\left( { - {I_1}} \right)}^N} + {I_1}} \right)}^2} - 4\ {\xi _{0{\rm{\ }}}}Ch} }}{{2\ {\xi _0}}}{\rm{\ }}. \end{eqnarray*}$$

(25c)

The positive radical sign provides correct solution for |$Ch\ = {\rm{\ }}0$|⁠. Figs B1 and B3 show sensitivity of the shape and size of the yield envelope to the change of the model parameters D₁, N and Ch values calculated using (25c). For N = 1 and non-cohesive material (Ch = 0) the envelope increases with decrease D₁ value keeping roughly self-similar shape (Fig. B1). Changing power index value N affects the shape of the envelope (Fig. B2). In these cases, D₁ values were rescaled to get the same volumetric strain 3.5% according to eq. ( 25b). Non-zero Ch value for the cohesive material with N = 1 and D₁ = 15 shifts envelope to larger differential strain values under low volumetric strains (Fig. B3). The shown set of surfaces for Ch = 0–1.0 × 10⁻⁶ corresponds to cohesive force values changing from zero to ∼15 MPa for rock with shear modulus |${\mu _0} = \ 10$| GPa. For typical cohesive force values of a few MPa, this model modification is important under relatively low confining pressures, but may be neglected when confining pressures are higher (of the order of tens of MPa).

Figure B1.

Yield envelope size increase with decrease in D-value (N = 1).

Figure B2.

Shape of the yield envelope for N-values changing between 0.5 and 2.0.

Figure B3.

The yield envelope is shifted up to higher differential strain values with increased Ch value (N = 1 and D = 15).

APPENDIX C

Anisotropic poroelastic damage model

The energy function of the anisotropic material with the damage tensor, |${{\rm{\Omega }}_{ij}}$|⁠, cannot be formulated only in terms of strain invariants I₁, I₂ of the elastic strain tensor, |${{\rm{\varepsilon }}_{ij}}$| ( |${I_1} = {{\rm{\varepsilon }}_{ii}}\ ,\ \ \ \ {I_2} = {{\rm{\varepsilon }}_{ij}}\ \ {{\rm{\varepsilon }}_{ij}}$|⁠) used for the scalar model (eq. 15). Following Murti et al. (1991) and Zhang & Cai (2010) the energy function should also depend on the invariants of the tensor, |$\varepsilon _{ij}^{( {\rm{\Omega }} )}$|⁠, which is the symmetrized product of the elastic strain and damage tensors:

$$\begin{eqnarray*} \varepsilon _{ij}^{\left( {\rm{\Omega }} \right)} = \frac{1}{2}\ \ \left( {{\varepsilon _{ik}}{{\rm{\Omega }}_{kj}} + {\varepsilon _{jk}}{{\rm{\Omega }}_{ki}}} \right). \end{eqnarray*}$$

(26)

The invariants |$I_1^{( {\rm{\Omega }} )}$| and |$I_2^{( {\rm{\Omega }} )}$| of this tensor are:

$$\begin{eqnarray*} I_1^{\left( {\rm{\Omega }} \right)} = {\varepsilon _{ik}}{{\rm{\Omega }}_{kj}}\ \ \ \ {\delta _{ij}} = {\varepsilon _{11}}{{\rm{\Omega }}_{11}}\ \ + {\varepsilon _{22}}{{\rm{\Omega }}_{22}} + {\varepsilon _{33}}{{\rm{\Omega }}_{33}} + 2{\varepsilon _{12}}{{\rm{\Omega }}_{12}} + 2{\varepsilon _{23}}{{\rm{\Omega }}_{23}} + 2{\varepsilon _{13}}{{\rm{\Omega }}_{13}} \end{eqnarray*}$$

(27)

$$\begin{eqnarray*} I_2^{\left( {\rm{\Omega }} \right)} &=& \varepsilon _{ij}^{\left( {\rm{\Omega }} \right)}\ \ \varepsilon _{ij}^{\left( {\rm{\Omega }} \right)}\ \ = {\left( {{\varepsilon _{11}}{{\rm{\Omega }}_{11}} + {\varepsilon _{12}}{{\rm{\Omega }}_{12}} + {\varepsilon _{13}}{{\rm{\Omega }}_{13}}} \right)^2} + {\left( {{\varepsilon _{12}}{{\rm{\Omega }}_{12}} + {\varepsilon _{22}}{{\rm{\Omega }}_{22}} + {\varepsilon _{23}}{{\rm{\Omega }}_{23}}} \right)^2} + {\left( {{\varepsilon _{13}}{{\rm{\Omega }}_{13}} + {\varepsilon _{23}}{{\rm{\Omega }}_{23}} + {\varepsilon _{33}}{{\rm{\Omega }}_{33}}} \right)^2}\nonumber \\ && + {\left( {{\varepsilon _{11}}{{\rm{\Omega }}_{12}} + {\varepsilon _{22}}{{\rm{\Omega }}_{12}} + {\varepsilon _{23}}{{\rm{\Omega }}_{13}} + {\varepsilon _{13}}{{\rm{\Omega }}_{23}} + {\varepsilon _{12}}\left( {{{\rm{\Omega }}_{11}} + {{\rm{\Omega }}_{22}}} \right)} \right)^2} + {\left( {{\varepsilon _{23}}{{\rm{\Omega }}_{12}} + {\varepsilon _{11}}{{\rm{\Omega }}_{13}} + {\varepsilon _{33}}{{\rm{\Omega }}_{13}} + {\varepsilon _{12}}{{\rm{\Omega }}_{23}} + {\varepsilon _{13}}\left( {{{\rm{\Omega }}_{11}} + {{\rm{\Omega }}_{33}}} \right)} \right)^2}\nonumber \\ && + {\left( {{\varepsilon _{13}}{{\rm{\Omega }}_{12}} + {\varepsilon _{12}}{{\rm{\Omega }}_{13}} + {\varepsilon _{22}}{{\rm{\Omega }}_{23}} + {\varepsilon _{33}}{{\rm{\Omega }}_{23}} + {\varepsilon _{23}}\left( {{{\rm{\Omega }}_{22}} + {{\rm{\Omega }}_{23}}} \right)} \right)^2} . \end{eqnarray*}$$

We extend energy function of the Panteleev et al. (2021) for the non-linear anisotropic damage rheology model by additional Biot terms similar to (15):

$$\begin{eqnarray*} F\ = \left[ {\frac{{{\lambda _0}}}{2}I_1^2 + {\mu _0}{I_2}} \right]{\rm{\ \ }} + \gamma \left[ {{\xi _0}I_2^{\left( {\rm{\Omega }} \right)} - I_1^{\left( {\rm{\Omega }} \right)}\sqrt {I_2^{\left( {\rm{\Omega }} \right)}} } \right] + \frac{1}{2}M{\left[ {\beta {I_1} - \zeta + {{\rm{\psi }}_{ij}}{\delta _{ij}}} \right]^2} + Ch{\rm{\ }}\frac{\gamma }{2}{\rm{\ }}{{\rm{\Omega }}_{ij}}{{\rm{\Omega }}_{ij}}. \end{eqnarray*}$$

(28)

The definitions (26)–(28) mean that the energy is the second order function of both elastic strain and damage tensors. Substituting energy function (28) into (5) and (6) the stress tensor is expressed as:

$$\begin{eqnarray*} {\sigma _{ij}} = \frac{{\partial F}}{{\partial {\varepsilon _{ij}}}}\ \ = {\lambda _0}\ \ {I_1}{\delta _{ij}} + 2{\mu _0}{\varepsilon _{ij}} - \gamma \frac{{\partial I_1^{\left( \Omega \right)}}}{{\partial {\varepsilon _{ij}}}}\sqrt {I_2^{\left( {\rm{\Omega }} \right)}} + \gamma \left( {{\xi _0} - \frac{{I_1^{\left( {\rm{\Omega }} \right)}}}{{2\sqrt {I_2^{\left( {\rm{\Omega }} \right)}} }}} \right)\frac{{\partial I_2^{\left( {\rm{\Omega }} \right)}}}{{\partial {\varepsilon _{ij}}}} + \beta M\ \left( {\beta {I_1} - \zeta + {{\rm{\psi }}_{kn}}{\delta _{kn}}} \right){\delta _{ij}}, \end{eqnarray*}$$

(29)

where

$$\begin{eqnarray*} \frac{{\partial I_1^{\left( \Omega \right)}}}{{\partial {\varepsilon _{lm}}}} = {{\rm{\Omega }}_{lm}}{\rm{\ }}, \end{eqnarray*}$$

$$\begin{eqnarray*} \frac{{\partial I_2^{\left( \Omega \right)}}}{{\partial {\varepsilon _{lm}}}} = \frac{1}{2}\ \ \left( {{\Omega _{li}}{\varepsilon _{ik}}{\Omega _{km}} + {\Omega _{mi}}{\varepsilon _{ik}}{\Omega _{kl}}} \right) + \frac{1}{2}\left( {{\Omega _{lj}}{\varepsilon _{mk}}{\Omega _{kj}} + {\Omega _{mj}}{\varepsilon _{lk}}{\Omega _{kj}}} \right). \end{eqnarray*}$$

The fluid pressure is almost identical to (17)

$$\begin{eqnarray*} {p_f} = \frac{{\partial F}}{{\partial \zeta }}\ \ = \ M\ \left( { - \beta {I_1} + \zeta - {{\rm{\psi }}_{kn}}{\delta _{kn}}} \right). \end{eqnarray*}$$

(30)

Similar to (18) the effective stress is defined as:

$$\begin{eqnarray*} {\sigma _{ij}} - \frac{{\partial F}}{{\partial {{\rm{\psi }}_{ij}}}}\ = \sigma _{ij}^{\rm eff}{\rm{\ }}. \end{eqnarray*}$$

(31)

In order to evaluate the energy dissipation and write the coupled damage-compaction kinetic relations (12), the damage-induced energy change should be calculated:

$$\begin{eqnarray*} \frac{{\partial F}}{{\partial {\Omega _{ij}}}} = \ \gamma \left[ {\left( {{\xi _0} - \frac{{I_1^{\left( {\rm{\Omega }} \right)}}}{{2\sqrt {I_2^{\left( {\rm{\Omega }} \right)}} }}} \right)\frac{{\partial I_2^{\left( {\rm{\Omega }} \right)}}}{{\partial {\Omega _{ij}}}} - {\varepsilon _{ij}}\sqrt {I_2^{\left( {\rm{\Omega }} \right)}} + Ch\ {{\rm{\Omega }}_{ij}}} \right]. \end{eqnarray*}$$

(32)

Formulating the tensor kinetic coefficients, we keep in mind that the coupling kinetic coefficient |${D_{ijkn}} = C_{ijnm}^{\Omega \psi }\ \ \ $| of eq. (13) should decrease with damage accumulation and expect that this will allow reproducing the directional Kaizer Effect observed under true triaxial loading (Browning et al. 2018). Similarly to the exponential relation (22), we suggest that the matrix |${D_{ijkn}}$| of the kinetic coefficients exponentially decrease with damage accumulation. Extending this idea to the complete tensor form and using the same dependency of the coupling term as in the scalar model (22), we suggest the following form of the coupling kinetic coefficient |${D_{ijkn}}$|⁠:

$$\begin{eqnarray*} {D_{ijkn}} = {D_1}\ \ \exp\left[ { - {D_2}{{\rm{\Omega }}_{ij}}} \right]\ {\delta _{kn}}\ {\left( { - {I_1}} \right)^N}\sqrt {{I_2}} , \end{eqnarray*}$$

(33)

where we use the standard definition of the exponent of the tensor X by means of its series representation (Hirsch et al. 1974)

$$\begin{eqnarray*} \exp\ \ \left( {\boldsymbol{X}} \right) = \mathop \sum \limits_{m = 0}^\infty \frac{{{{\boldsymbol{X}}^m}}}{{m!}}{\rm{\ }}. \end{eqnarray*}$$

(34)

Note that the principal values of the tensor |$\exp( {\boldsymbol{X}} )$| are equal to the |$\exp( {{X_k}} )$|⁠.

The most conservative assumption to define the components of the matrix |$C_{ijnm}^{\Omega \Omega }$| in (13) is the absence of the interaction between different components, that is |$C_{ijnm}^{\Omega \Omega } = \ L\ ( {{\delta _{ik}}\ {\delta _{jn}} + {\delta _{in}}\ {\delta _{jk}}} )$|⁠, and L is proportional to |$1/\sqrt {{{\rm{\Omega }}_{ij}}{{\rm{\Omega }}_{ij}}} $|⁠, as it was suggested by Panteleev et al. (2021) using results of true-triaxial rock mechanics experiments ignoring effects of compaction. With these kinetic coefficients, the eq. (12b) for the damage evolution has the same structure as the damage kinetic eq. (23) of the scalar model.

$$\begin{eqnarray*} \frac{{{\rm d}{\Omega _{ij}}}}{{{\rm d}t}} = {D_1}\ \ \exp\left[ { - {D_2}{{\rm{\Omega }}_{ij}}} \right]{\rm{\ \ }}{\left( { - {I_1}} \right)^N}\sqrt {{I_2}} {\rm{\ }}\sigma _m^{\rm eff} + {\rm{\ }}\frac{L}{{\sqrt {{{\rm{\Omega }}_{kn}}{{\rm{\Omega }}_{kn}}} }}\left[ {{\varepsilon _{ij}}\sqrt {I_2^{\left( {\rm{\Omega }} \right)}} + \left( {\frac{{I_1^{\left( {\rm{\Omega }} \right)}}}{{2\sqrt {I_2^{\left( {\rm{\Omega }} \right)}} }} - {\xi _0}} \right)\frac{{\partial I_2^{\left( {\rm{\Omega }} \right)}}}{{\partial {\Omega _{ij}}}} - Ch{\rm{\ }}{{\rm{\Omega }}_{ij}}} \right]. \end{eqnarray*}$$

(35)

The equation for the evolving compaction (12a) includes two terms. The first term, equal to |$C_{ijnm}^{\psi \psi }\sigma _{nm}^{\rm eff}$| describe compaction/extension with the rate proportional to the effective stress. The second one is the damage-related coupling term, which is the |$\frac{{\partial F}}{{\partial {\Omega _{ij}}}}$| multiplied to |${D_{ijkn}}$|⁠, and describes the compaction or dilation induced by the damage growth. Representing the effective stress as a superposition of the volumetric (⁠|${P_e}$|⁠) and deviatoric (⁠|${\tau _{ij}}$|⁠) components allows to describe different mechanisms of the irreversible strain accumulation or |${{\rm{\psi }}_{ij}}$| kinetics. The pressure driven compaction under hydrostatic load is described by a well-known porosity reduction to its pressure-dependent equilibrium value or Athy's (1930) law. Lyakhovsky et al. (2022) modified the scalar Athy relation accounting not only for volumetric effects, but also directional effects associated with changes of the shape of pore space and accumulation of irreversible deviatoric strain components. The suggested 3-D equilibrium compaction strain, |${\rm{\psi }}_{ij}^{( {\rm eq} )}$|⁠, depends on both pressure and deviatoric stress components:

$$\begin{eqnarray*} {\rm{\psi }}_{ij}^{\left( {\rm eq} \right)} = {B_0}\ \ \left[ {{\delta _{ij}} - \ \exp\left( { - \frac{{\sigma _m^{\rm eff}}}{{{B_1}}}{\delta _{ij}} - \frac{{{\tau _{ij}}}}{{{B_2}}}} \right)} \right] \end{eqnarray*}$$

(36)

and suggested kinetics of the pressure-driven 3-D compaction in the form

$$\begin{eqnarray*} \frac{{{\rm d}{{\rm{\psi }}_{ij}}}}{{{\rm d}t}} = \ A\ \left( {{\rm{\psi }}_{ij}^{\left( {\rm eq} \right)} - {{\rm{\psi }}_{ij}}} \right)\ \sigma _m^{\rm eff}. \end{eqnarray*}$$

(37)

Neglecting the term with deviatoric stress or taking |${B_2} \to \infty $| in (36) reduces both equilibrium compaction (36) and kinetic eq. (37) to the traditional scalar Athy's (1930) law formulated in terms of material porosity.

Experimental studies suggest that permanent inelastic deformation is accumulated in high porosity rocks, but also is related to the damage accumulation. It starts accumulating with the onset of the AE and increases all the way up to brittle failure (e.g. Lockner 1993; Martin & Chandler 1994, 1998). This process is usually associated with the growth of microcracks and frictional sliding between grains, rather than closure of voids or opening space between grans. For similar reasons, Hamiel et al. (2004a) related the rate of irreversible strain accumulation with the rate of their scalar damage growth. Keeping in mind that the scalar damage variable is equivalent to the squared damage tensor, we extend their relation to the tensor form keeping the same structure of functional relations:

$$\begin{eqnarray*} \frac{d}{{{\rm d}t}}\ \ \varepsilon _{ij}^{\left( {ir} \right)} = \left \{ \begin{array}{@{}*{1}{c}@{}} {{C_V}\frac{{d\left( {{{\rm{\Omega }}_{in}}{{\rm{\Omega }}_{jm}}} \right)}}{{{\rm d}t}}{\tau _{nm}},\ \ \ \ \frac{{d\left( {{{\rm{\Omega }}_{in}}{{\rm{\Omega }}_{jm}}} \right)}}{{{\rm d}t}} \gt 0}\\ {0,\ \ \ \ \ \ \ \ \ \ \ \frac{{d\left( {{{\rm{\Omega }}_{in}}{{\rm{\Omega }}_{jm}}} \right)}}{{{\rm d}t}} \le 0}. \end{array}\right. \end{eqnarray*}$$

(38)

Here we combine both mechanisms and represent the static value of the tensor |${\rm{\psi }}_{ij}^{( S )}$| as a sum of

$$\begin{eqnarray*} {\rm{\psi }}_{ij}^{\left( S \right)} = {\rm{\psi }}_{ij}^{\left( {\rm eq} \right)}\ \ + \varepsilon _{ij}^{\left( {ir} \right)}. \end{eqnarray*}$$

(39)

The final kinetic eq. (12a) for the |${\psi _{ij}}$| tensor incorporates all the discussed mechanisms (to avoid very length equation, |${D_{ijkn}}$| and |$\frac{{\partial F}}{{\partial {\Omega _{kn}}}}$| are not substituted here):

$$\begin{eqnarray*} \frac{{{\rm d}{\psi _{ij}}}}{{{\rm d}t}} = \ A\ \left( {{\rm{\psi }}_{ij}^{\left( {\rm eq} \right)} - {{\rm{\psi }}_{ij}}} \right)\ \sigma _m^{\rm eff} + {C_V}\frac{{d\left( {{{\rm{\Omega }}_{in}}{{\rm{\Omega }}_{jm}}} \right)}}{{{\rm d}t}}{\tau _{nm}} - {D_{ijkn}}\frac{{\partial F}}{{\partial {\Omega _{kn}}}}. \end{eqnarray*}$$

(40)

The kinetic equation for damage (35) and compaction (40) provide the close system of equations defining the 3-D evolution of the material properties. To study some basic model properties, we follow the assumption of Panteleev et al. (2021) that the principal directions of the damage tensor match the orientation of the principal loading axes. This assumption is supported by results of true-triaxial experiments (Browning et al. 2017, 2018) demonstrating orientation of distributed micro-cracks in the sample volume under limited load, well before their localization into narrow fault zones. Adopting the above assumptions, we rewrite the stress–strain and damage kinetic relations for principal values in k-direction (no summation):

$$\begin{eqnarray*} {\sigma _k} = {\lambda _0}\ \ \ {I_1} + 2{\mu _0}{\varepsilon _k} + 2{\mu _1}{\rm{\Omega }}_k^2{\varepsilon _k} - \gamma {{\rm{\Omega }}_k}\sqrt {I_2^{\left( {\rm{\Omega }} \right)}} - \gamma {\xi ^{\left( {\rm{\Omega }} \right)}}{\rm{\Omega }}_k^2{\varepsilon _k} + \beta \ {p_f} \end{eqnarray*}$$

(41)

and

$$\begin{eqnarray*} \frac{{{\rm d}{\Omega _k}}}{{{\rm d}t}} = {D_1}\ \ \exp\left[ { - {D_2}{{\rm{\Omega }}_k}} \right]{\rm{\ }}{\left( { - {I_1}} \right)^N}\sqrt {{I_2}} {\rm{\ }}\sigma _m^{\rm eff} + \frac{L}{{\sqrt {{{\rm{\Omega }}_l}{{\rm{\Omega }}_l}} }}\left[ {{\varepsilon _k}\sqrt {I_2^{\left( {\rm{\Omega }} \right)}} + \left( {\frac{{I_1^{\left( {\rm{\Omega }} \right)}}}{{\sqrt {I_2^{\left( {\rm{\Omega }} \right)}} }} - 2{\xi _0}} \right)\varepsilon _k^2{{\rm{\Omega }}_k} - Ch{\rm{\ }}{{\rm{\Omega }}_k}} \right]. \end{eqnarray*}$$

(42)

We note that for the isotropic damage ( |${{\rm{\Omega }}_1} = {{\rm{\Omega }}_2}\ \ = {{\rm{\Omega }}_3}\ $|⁠) summation of free damage rate components give the same equation as for the scalar damage with the same damage-dependent yield envelope.

$$\begin{eqnarray*} \frac{{{\rm d}{\psi _k}}}{{{\rm d}t}} = \ A\ \left( {{\rm{\psi }}_k^{\left( {\rm eq} \right)} - {{\rm{\psi }}_k}} \right)\ \sigma _m^{\rm eff} + {C_V}\frac{{d\left( {\Omega _k^2} \right)}}{{{\rm d}t}}{\tau _k} + {D_1}\exp\left[ { - {D_2}{{\rm{\Omega }}_k}} \right]\ {\left( { - {I_1}} \right)^N}\sqrt {{I_2}} \ \left[ {{\varepsilon _k}\sqrt {I_2^{\left( {\rm{\Omega }} \right)}} + \left( {\frac{{I_1^{\left( {\rm{\Omega }} \right)}}}{{\sqrt {I_2^{\left( {\rm{\Omega }} \right)}} }} - 2{\xi _0}} \right)\varepsilon _k^2{{\rm{\Omega }}_k} - Ch{\rm{\ }}{{\rm{\Omega }}_k}} \right]. \end{eqnarray*}$$

(43)

The first term of eq. (43) represents the compaction prior to the onset of the damage accumulation, which is 3-D extension of the scalar Athy's law compaction. According to this term the compaction approaches to its stress-dependent equilibrium value with the rate proportional to the effective pressure. The second term is 3-D equivalent to the damage-dependent irreversible strain accumulation with inverse of the effective viscosity or fluidity proportional to the rate of the damage accumulation. This term describes extension or compaction depending on the sign of the deviatoric stress component. The last term represents the coupling between damage and porosity kinetics. Its sign, extension or compaction, is defined by the expression in the square brackets and depends on the loading and damage values.