Sliding Crack Model for Nonlinearity and Hysteresis in the Triaxial Stress‐Strain Curve of Rock, and Application to Antigorite Deformation

Under triaxial deviatoric loading at stresses below failure, rocks generally exhibit nonlinearity and hysteresis in the stress‐strain curve. In 1965, Walsh first explained this behavior in terms of frictional sliding along the faces of closed microcracks. The hypothesis is that crack sliding is the dominant mode of rock inelasticity at moderate compressive stresses for certain rock types. Here we extend the model of David et al. (2012, https://doi.org/10.1016/j.ijrmms.2012.02.001) to include (i) the effect of the confining stress; (ii) multiple load‐unload cycles; (iii) calculation of the dissipated strain energy upon unload‐reload; (iv) either frictional or cohesive behavior; and (v) either aligned or randomly oriented cracks. Closed‐form expressions are obtained for the effective Young's modulus during loading, unloading, and reloading, as functions of the mineral's Young's modulus, the crack density, the crack friction coefficient and cohesion for the frictional and cohesive sliding models, respectively, and the crack orientation in the case of aligned cracks. The dissipated energy per cycle is quadratic and linear in stress for the frictional and cohesive models, respectively. Both models provide a good fit to a cyclic loading data set on polycrystalline antigorite, based on a compilation of literature and newly acquired data, at various pressures and temperatures. At high pressure, with increasing temperature, the model results reveal a decrease in friction coefficient and a transition from a frictionally to a cohesively controlled behavior. New measurements of fracture toughness and tensile strength provide quantitative support that inelastic behavior in antigorite is predominantly caused by shear crack sliding and propagation without dilatancy.


Introduction
It is well known that the mechanical behavior of polycrystalline brittle rocks under confined compressive loading is to a great extent controlled by the presence of crack-like flaws or voids located within grains and along grain boundaries. This applies to both processes of elastic and inelastic deformation, as well as rock failure (Paterson & Wong, 2005). These features are also observed in other brittle materials such as ceramic composites (Marshall & Oliver, 1987) or concrete (Shah et al., 1995).
We examine here the hypothesis that, under deviatoric loading, inelastic deformation is predominantly accommodated by shear-induced sliding of preexisting microcracks, without in-or out-of-plane crack growth. This hypothesis seems valid for stress conditions and rock types as follows: (a) The confining pressure should be sufficiently high so that most preexisting microcracks are closed (typically a few hundred of MPa in rocks) with their surfaces in contact. (b) The compressive loading stress should be below that required for the onset of crack propagation. (c) Rock types are low-porosity rocks that contain microcracks or, more generally, "planar" surfaces amenable for sliding under shear. (d) Crack propagation and notably the possible onset of dilatancy seem to occur at stresses substantially greater than the yield point.
The recent cyclic loading experiments of David et al. (2018) on axially loaded polycrystalline antigorite provide motivation for this hypothesis. To about 90% of the failure stress, the load-unload stress-strain curve (e.g., Figure 4e of that publication) can be divided into four regimes and interpreted as follows. Initially, during loading at low stress, the behavior is linear elastic with a Young's modulus equal to the "intrinsic" modulus, that is, that of the uncracked solid. At higher stress, the behavior deviates from linear elasticity and becomes nonlinear, with a Young's modulus that is stress dependent and reduced compared to that of the uncracked solid, owing to sliding on the crack interfaces (without crack growth). At the beginning of unloading, the behavior is linear with a Young's modulus again equal to the intrinsic modulus, which reflects a 10.1029/2019JB018970 transient delay in crack "reverse" or "back" sliding, producing hysteresis in the stress-strain curve. At lower stresses during unloading, the onset of reverse sliding on cracks is associated with nonlinear behavior, with a tangent modulus that is again stress dependent and lower than the intrinsic modulus. The initial delay in activating reverse sliding on the cracks causes the tangent modulus to be even lower than during forward sliding. However, the strain does not return to the origin when stress is removed, and a non-negligible amount of "permanent" strain is created. The measurements of volumetric strain and seismic velocities by David et al. (2018), combined with microstructural observations, reveal that such behavior is observed without detectable dilatancy or crack growth. In addition, such observations are made at the highest available confining pressure (150 MPa) which appears to be sufficient to close most microcracks (David et al., 2018). To a broader scope, the hypothesis of "Mode II," nondilatant shear crack sliding seems relevant to the inelastic behavior of phyllosilicates generally, such as lizardite (Escartín et al., 1997) or talc (Escartín et al., 2008). The transition from "Mode I" to "Mode II" microcracking at increasing pressure is also considered as an important mechanism involved in the brittle-plastic transition, for example, as observed in quartz (Hirth & Tullis, 1994). The first analytical model quantitatively accounting for the effect of frictional sliding on closed microcracks in rocks was proposed by McClintock and Walsh (1962) to modify Griffith's theory for rock failure, followed by a fundamental paper by Walsh (1965) who considered both crack closure and frictional sliding on closed cracks, and calculated the stress-dependent Young's modulus and stress-strain curve during unaxial loading for a rock containing randomly oriented cracks. Useful analysis of backsliding on cracks during unloading was also given in that paper. Kachanov (1982a) proposed a rigorous three-dimensional analysis on a rock containing closed microcracks and extended Walsh's model to triaxial state of stress. However, no quantitative analysis of the rock behavior during unloading was given. Horii and Nemat-Nasser (1983) considered both crack closure and crack frictional sliding under triaxial stress, this time accounting for crack interactions by using the self-consistent theory, but again did not consider the unloading process. Lehner and Kachanov (1995) extended Kachanov's analysis to the unloading and reloading process, but their model was not tested against experimental data. A complete analysis of unloading and reloading of a material containing closed microcracks was given by Lawn and Marshall (1998). They added a cohesive term to the frictional constitutive law for crack sliding. However, their solutions were only given for uniaxial compression, and no experimental data were analyzed. David et al. (2012) extended Walsh's model to an entire load-unload cycle during uniaxial compression by accounting for both crack closure and frictional sliding and fitted the model to stress-strain data on sandstone and granite. Bruno and Kachanov (2013) incorporated the effect of the nonclosable pores in addition to crack closure and frictional sliding during a uniaxial load-unload cycle. Their model, which uses approximate schemes to account for crack and pore interactions, was applied to stress-strain as well as microstructural data on ceramic. In parallel, a large number of micromechanical crack-based models for the inelasticity of brittle materials have coupled the "sliding crack mechanism" with either the dilatancy associated with the growth of Mode I tension cracks or kinks at the crack tips (Ashby & Sammis, 1990;Basista & Gross, 1998;Kachanov, 1982b;Lehner & Kachanov, 1996;Moss & Gupta, 1982;Nemat-Nasser & Horii, 1982;Nemat-Nasser & Obata, 1988;Stevens & Holcomb, 1980) or self-similar crack growth (Fanella & Krajcinovic, 1988;Gambarotta & Lagomarsino, 1993), some of which were successfully fit to experimental data. While the main emphasis of these models was on stress-induced crack growth, it is difficult to extract closed-form expressions for the simplified case of sliding without crack propagation from these publications. In addition, the case of the reloading of sliding cracks was not considered.
Hence, here we propose an analytical model for rock inelasticity under deviatoric loading, based purely on sliding of preexisting microcracks whose surfaces are initially in contact, and without crack propagation. The model is a direct extension of the model of David et al. (2012) to account for the effect of a constant lateral or "confining" stress. For generality and concision of analysis, we use a constitutive law for crack sliding that combines both a cohesive and a frictional "Coulomb-type" resistance to sliding. Following Lawn and Marshall (1998), incorporation of the effect of a confining stress is indeed identical to adding a cohesive term in the mechanics of the problem; the analytical results for the sliding crack model under triaxial stress can thus be derived from Lawn and Marshall (1998) using the superposition principle, as we shall see below. Two separate cases are considered: aligned or randomly oriented cracks. The emphasis of the model results and their application to fitting experimental data are then kept separate for two "end-member" models, purely frictional and purely cohesive sliding, for reasons related to simplicity of physical interpretation. For example, (i) in contrast to the case of uniaxial compression, if a confining stress is applied, there is no need to include a cohesion term in the constitutive law for crack sliding in order to predict yield-type behavior; (ii) we wish to minimize the number of parameters used to fit the experimental data; and (iii) "friction" and "cohesion" have different physical origins, and a cohesion-only model may be used to simulate intracrystalline slip, whereas friction (i.e., normal stress-dependent slip) should require a minimum degree of crack roughness and the presence of interface asperities at the microscale or nanoscale (Bowden & Tabor, 2001). In addition to giving closed-form expressions for the loading, unloading, and reloading portions of the stress-strain curve, we calculate the dissipated strain energy between unloading and subsequent reloading, for the purely frictional and purely cohesive sliding cases. To render the problem amenable to analytical treatment, we assume that the rock is initially isotropic and that crack interactions can be neglected. The latter assumption underestimates the effect of cracks on rock inelasticity (Horii & Nemat-Nasser, 1983) but seems quantitatively valid up to moderate cracks concentrations (e.g., Kachanov, 1982a). The model is then fit to cyclic loading stress-strain data on polycrystalline antigorite from David et al. (2018) at 150 MPa confining pressure and room temperature, and from newly conducted experiments at 1 GPa at room temperature, 400 • C and 500 • C in a Griggs apparatus.

Effect of a Single Crack During Cyclic Loading: Crack Displacement-Stress Relations
Consider a rock specimen of cross-sectional area A, containing a single crack of length 2c, that is closed so that the two opposite faces are in contact and can slide past each other (Figure 1; see nomenclature given in Table 1). The crack is subjected to a vertical stress 1 and a lateral stress 2 . The convention used is that compressive stresses are positive. The "confining" or minimum stress 2 is held constant, while under deviatoric loading the maximum compressive stress 1 ≥ 2 . The normal to the crack plane makes an angle to 1 . The resolved normal and shear stresses on the crack are then related to the stresses as n = 1 cos 2 + 2 sin 2 and = ( 1 − 2 ) cos sin , respectively. In this section we determine the average displacement on the crack, b k , to the applied state of stress, for the cases k = L (loading), k = U (unloading), and k = RL (reloading). A fundamental relation is that the sliding displacement increases linearly with a driving "effective sliding stress," k eff (where k = L; U; RL) as (Stevenson, 1945) where E 0 is the Young's modulus of the uncracked rock, and where the effective sliding stress k eff is defined as the difference between an applied shear stress driving sliding and a shear stress resisting sliding. k eff must be positive for a crack to slide. Its expression depends on the crack orientation, the acting stresses, and the constitutive parameters for sliding, but also on the loading history (loading, unloading, reloading) as detailed below. Note that Equation 1 is written under the assumption of plane stress which is the one used here. Results can be transformed into one appropriate for plane strain by replacing E 0 by E 0 ∕(1 − 2 0 ), where 0 is the Poisson's ratio of the uncracked rock. In any event, the term (1 − 2 0 ) is very close to unity.

Loading
Consider the model of a sliding crack endowed with both cohesive and frictional "Coulomb-type" resistance. Sliding is resisted by a frictional shear stress, f , expressed as (Nemat-Nasser & Obata, 1988) f = y + n = y + ( 1 cos 2 + 2 sin 2 ), where y is a constant "cohesive" yield stress and denotes the friction coefficient. A simple substitution is needed to make the analysis of Lawn and Marshall (1998) applicable for deviatoric loading with vertical This equation is in the form of a constant ( y + 2 ) plus Coulomb friction and is thus identical to the constitutive law in Lawn and Marshall (1998) for uniaxial loading ( f = c + sin 2 , where is the uniaxial stress), with their cohesive stress c replaced by y + 2 and a different angle definition, = ∕2 − . Therefore, with the substitutions c → y + 2 and → 1 − 2 , the analysis of Lawn and Marshall (1998) can be directly applied to deviatoric loading.
The effective sliding stress on crack during loading (L), L eff , is the difference between the resolved shear stress and the frictional resistive shear stress: The condition for crack sliding is that L eff > 0. This occurs when 1 > L y , where the "crack yield stress during loading," L y , is found by solving for the stress at which L eff = 0 and expressed as L y = y + 2 (cos sin + sin 2 ) cos sin − cos 2 .

10.1029/2019JB018970
We note that the model does not necessarily require a cohesion term in the constitutive law to predict yield-type behavior; here such behavior already arises simply from the application of a confining stress.
Once sliding is activated, by differentiating the crack displacement-stress Equation 1 and the expression (4) for the effective stress during loading, an increment of crack sliding displacement db L is related to a stress increment d 1 as hence, the crack displacement during loading is given by

Unloading
Now consider that the rock has been loaded to a maximum stress ⋆ 1 (where the superscript ⋆ generally denotes the value of a variable taken at the maximum stress) and that a sliding displacement has been accumulated on the crack. During unloading, the restoring force for backsliding is the elastic strain accumulated at the crack tips during loading. The effective stress driving backsliding is the difference between the effective stress at the maximum stress, ⋆ eff = L eff ( 1 = ⋆ 1 ), which provides the restoring force, and the joint action of the frictional resisting stress and the resolved applied shear stress, which both act against backsliding (Nemat-Nasser & Obata, 1988). Using Equation 4 and previous definitions, it is found that The condition for backsliding is U eff > 0. At the onset of unloading ( 1 = ⋆ 1 ), when the direction of loading is reversed, using the definition of U eff above we obtain that U eff = −2 ⋆ f < 0. Hence, regardless of its orientation, a crack that has been sliding remains "stuck" at the beginning of unloading. The load must decrease by a finite amount before backsliding is activated. Backsliding occurs when 1 < U y , where the "crack yield stress during unloading" U y is found by solving for the stress at which U eff = 0 and is given by U y = ⋆ 1 (cos sin − cos 2 ) − 2( y + 2 sin 2 ) cos sin + cos 2 .
Similarly to the loading case, by differentiating Equations 1 and 8, and under the algebraic convention that a positive displacement increment means forward sliding, once backsliding is activated an increment of crack backsliding displacement db U is related to a stress increment as If b ⋆ denotes the maximum displacement on crack at the end of loading (i.e., from Equation 7, b ⋆ = ( c∕E 0 )(cos sin − cos 2 )( ⋆ 1 − L y )), the crack displacement during unloading is then given by After some algebra, by evaluating b U at 1 = 2 in Equation 11 and by using Equations 5 and 9, we find that the residual displacement on the crack at the end of unloading is equal to ( c/E 0 )( y + 2 ). Hence, the residual displacement on a crack at the end of a load-unload cycle is not only independent of the crack orientation but also on the magnitude of the applied stress, a remarkable result that was not a priori intuitive.

Reloading
Assume now that the rock has been loaded to a maximum stress ⋆ 1 and unloaded completely so that 1 = 2 . In this section we demonstrate that the behavior of the crack during a subsequent reloading depends on its accumulated history during the load-unload cycle, and how. The residual shear stress and displacement on the crack at the end of the load-unload cycle must therefore be calculated. This allows to obtain the effective sliding stress on the crack during reloading, RL eff , the yield stress during reloading, RL y , and the cracked displacement during reloading, b RL . Based on previous analysis, at the end of unloading the crack can be in three possible configurations: (i) No sliding during the previous load-unload cycle. This occurs if ⋆ 1 < L y . The shear stresses are all relaxed on the crack face, so in this case RL eff = L eff , and RL y = L y . The crack displacement-stress relation during reloading is given by Equation 7. (ii) Forward sliding only during the previous load-unload cycle. The joint condition for this to occur is that ⋆ 1 > L y (sliding) and U y < 2 (no backsliding). The shear stress accumulated on the crack at the end of loading, ⋆ eff , is not relaxed by any backsliding during unloading and therefore must be overcome during reloading for the crack to slide again the forward direction. The effective shear stress driving sliding during reloading, in this situation, is thus expressed as Therefore, the crack yield stress during reloading, found by solving for the stress at which RL eff = 0, is simply RL y = ⋆ 1 . Hence, regardless of the crack orientation, a crack that had been sliding, but not did not experience backsliding during a load-unload cycle remains "stuck" until the stress during reloading exceeds the maximum stress during the previous cycle. By differentiating Equations 1 and 12, we find that once sliding resumes an increment of crack sliding displacement during reloading, db RL , is related to a stress increment as in Equation 6 for loading. Hence, according to previous analysis, the crack displacement during reloading is simply equal to b ⋆ when is 1 < ⋆ 1 and described by Equation 7 when 1 ≥ ⋆ 1 . (iii) Backsliding during the previous load-unload cycle. The sufficient condition for this to occur is that U y > 2 . In this case the residual shear stress on the crack at the end of unloading is equal to the difference between the shear stress accumulated during loading (i.e., ⋆ eff ) and the amount of shear stress relaxed by backsliding at the end of unloading (i.e., U eff ( 1 = 2 )). Using Equations 4 and 8, respectively, such a difference is equal to y + 2 . As previously found for the residual displacement, the residual (or unrelaxed) shear stress on crack at the end of a load-unload cycle is independent of its orientation, but also on the maximum applied stress. Since this residual stress must be overcome during reloading for the crack to slide again in the forward direction, the effective shear stress driving sliding during reloading will be given by The crack yield stress during reloading, found by solving for the stress at which RL eff = 0, is RL y = 2 y + 2 [cos sin + (1 + sin 2 )] cos sin − cos 2 .
It is worth comparing the crack yield stress during reloading to that during loading, but also to the maximum stress during the previous load-unload cycle. First, the yield stress on crack during reloading (14) is higher than that for loading (5), as could be expected. Additional comparison of the analytical expressions for the yield stress during unloading (9) and that for reloading (14) reveals, after some algebra, that the condition for backsliding ( U y > 2 ) is entirely equivalent to RL y < ⋆ 1 .
Hence, it is always the case that L y < RL y < ⋆ 1 . Note that the condition for the occurrence of backsliding conveniently written as (15) is used thereafter. A crack having undergone backsliding during a load-unload Figure 2. Crack displacement-stress diagrams for a body in compression containing a single closed shear crack, for the purely frictional (a-c) and purely cohesive sliding (d-f) cases. Solid lines: loading segments; dashed lines: unloading segments. Squares: yield points. In (c) and (f), the displacement-stress diagram is that of a body subjected to six successive load-unload cycles (OABC and CBBC to ⋆ 1 = 4 2 , CBDEF and FGDEF to ⋆ 1 = 7 2 , and FGHIF and FGHIF to ⋆ 1 = 10 2 ). cycle will remain "stuck" during reloading until 1 = RL y , but, contrary to the previous case (ii), sliding will resume on the crack at a stress below the maximum stress of the previous load-unload cycle. Finally, by differentiating Equations 1 and 13, we find that once sliding resumes an increment of crack sliding displacement during reloading is related to a stress increment as in Equation 6 for loading. Hence, according to previous analysis, the crack displacement during reloading is simply equal to ( c/E 0 )( y + 2 ) when 1 < RL y and described by Equation 7 when 1 ≥ RL y .

Representative Crack Displacement-Stress Diagrams for the Purely Frictional and Cohesive Sliding Models
The analytical results have been given above for a combined analysis of a crack endowed with both cohesive and frictional "Coulomb-type" resistance. Results for the purely frictional and the purely cohesive cases can be simply obtained by setting y = 0 and = 0, respectively, in the equations above. Figure 2 shows representative crack displacement versus stress diagrams for the purely frictional (a-c) and purely cohesive (d-f) sliding models. The applied stress 1 is normalized to the confining stress 2 . With this definition of normalized stress and the form taken by the crack displacement-stress relations (e.g., (7)), the crack displacement is normalized to c 2 /E 0 .
In the purely frictional case ( y = 0), Figure 2a shows the dependence of the crack displacement-stress response for three crack orientations ( = 38 • , 50 • , and 60 • ) and one fixed value of the friction coefficient ( = 0.6) for a body subjected to one load-unload cycle to ⋆ 1 = 10 2 . During both loading and unloading, the yield stress and also the crack displacement-stress rate are orientation dependent. In the case of a crack inclined at = 38 • , no backsliding occurs. For the other two cases ( = 50 • and = 60 • ) the crack displacement-stress rate during unloading is greater than for loading; however, the residual displacement left at the end of unloading is independent of the crack orientation. Figure 2b shows the equivalent plot for three values of the friction coefficient ( = 0.3, 0.6, and 0.9) and one fixed crack orientation ( = 60 • ), during one load-unload cycle to ⋆ 1 = 10 2 . A lower value of facilitates sliding in that sliding occurs at a lower stress and at a greater rate. When = 0.9, no backsliding occurs. For the two other cases ( = 0.6 and = 0.3), the residual displacement increases with . Figure 2c shows the crack displacement-stress 10.1029/2019JB018970 response for a body subjected to six successive load-unload cycles, for a given crack orientation ( = 60 • ) and friction coefficient ( = 0.6): two cycles to a maximum stress ⋆ 1 = 4 2 , followed by two cycles to ⋆ 1 = 7 2 , and by two cycles at ⋆ 1 = 10 2 . This illustrates the influence of increasing the loading stress and also the effect of a repeated load-unload cycle to the same stress. During the first cycle to ⋆ 1 = 4 2 (OABC), sliding occurs from the yield point (A) to the maximum loading stress (B); there is no backsliding during unloading (BC), and a residual displacement is left (at C). At the beginning of the second cycle to ⋆ 1 = 4 2 (CBBC), the reloaded crack is in the configuration (ii) described in section 2.3. Reactivation of sliding would thus require increasing the stress above the maximum stress of the previous cycle; hence, the crack remains "stuck" during both reloading (CB) and unloading (BC). The following cycle (CBDEF) is performed to a greater maximum stress than previously ( ⋆ 1 = 7 2 ). During reloading (CBD), sliding accordingly resumes at yield point B until the maximum loading stress is reached (D). During unloading (DEF), the crack is initially "stuck" (DE), after which backsliding occurs at the unloading yield point E. The residual displacement left (at F) is greater than during the previous cycle (at C). At the beginning of the second cycle to ⋆ 1 = 7 2 (FGDEF), the reloaded crack is this time in the configuration (iii) described in section 2.3. During reloading (FGD), reactivation of sliding occurs at the "reloading yield point" G, which is above that for initial sliding (A) but below the maximum stress of the previous cycle (D), after which the behavior is identical to that of previous cycle along the GDEF segment. The following cycle (FGHIF) is performed to a greater maximum stress than previously ( ⋆ 1 = 10 2 ). As the reloaded crack is initially in the configuration (iii), sliding resumes again at yield point G until the maximum loading stress is reached (H). During unloading (HIF), the crack is initially "stuck" (HI), after which backsliding occurs at the unloading yield point (I), which is above that of the previous cycle to a lower stress (E). However, the end of the unloading segment (EF) and hence the residual displacement is the same that during previous cycle. During the second cycle to ⋆ 1 = 10 2 (FGHIF), the behavior is entirely identical to the previous cycle.
The choice of the representative values of the parameters in the purely cohesive case ( = 0; Figures 2d-2f) is made to highlight the strong similarities of the crack behavior with that described in detail above for the frictional case. The main difference between the two models is that the slope of the crack displacement-stress for the cohesive sliding model is the same during loading and unloading and only depends on the crack orientation (Figures 2d and 2e). This is expected as once the cohesive yield stress y is overcome by the shear stress acting on the crack face, in the constitutive law for sliding there is no reason for the crack displacement-stress path to depend on y .
Overall, under triaxial loading, both purely frictional and cohesive sliding models predict yield-type behavior. Another main feature of the behavior is that, although the crack displacement upon reloading depends on the previous load-unload history, as soon as the stress reaches the maximum stress of the previous cycle ( ⋆ 1 ), the structure of the rock is entirely "reset." Hence, the rock behavior is reversible upon reloading; that is, the crack displacement at ⋆ 1 during reloading will be exactly the same as the one at ⋆ 1 during a previous cycle. Finally, once backsliding is activated, the residual displacement on a sliding crack does not depend on its orientation, nor the magnitude of the maximum stress during a load-unload cycle.

Body With Aligned Cracks 3.1.1. Loading
Consider first a single crack oriented at an angle to 1 that is sliding under stress during loading (i.e., 1 > L y ). In a representative area A, an increment of sliding displacement db L on the crack produces an increment of "inelastic" strain d i 1 in the direction of 1 as given by (Nemat-Nasser & Obata, 1988) Note that the lateral strain (in the direction of 2 ) can be simply calculated by replacing the term sin(2 ) by − cos(2 ) in this equation. Inserting the incremental crack displacement-stress relation (6) in previous equation, Now that the inelastic strain increment has been explicitly related to the applied stress, the total strain increment in the rock, d 1 , is the sum of the elastic strain and the inelastic strain increments: If we denote by E L the effective Young's modulus of the rock during loading, since by definition d 1 = d 1 /E L , inserting Equation 17 in Equation 18 we obtain, by eliminating d 1 , Now imagine that the rock contains an array of N cracks per representative area A, having the same length 2c and all oriented at the same angle to 1 . If we invoke the "no-interaction" approximation, the excess strain due to each sliding crack is computed as if it were an isolated crack; therefore, the crack contributions are purely additive. The c 2 /A term on the right-hand side of the previous equation then becomes N c 2 /A, and the effective Young's modulus during loading will be where L y is given by Equation 5 and Γ = Nc 2 ∕A is the two-dimensional crack density. Crack sliding causes a modulus deficit relative to the uncracked solid which is easily estimated from Equation 20. In the purely frictional case ( y = 0), typically, for a friction coefficient = 0.5 and a crack oriented at 45 • to 1 , E 0 /E L ≈ 1 + 0.8Γ. Taking a crack density Γ = 0.5, E L ≈ 0.7E 0 ; that is, the frictional sliding on the array of cracks causes a 30% modulus deficit relative to the uncracked solid. In the purely cohesive case ( = 0), Equation 20 becomes E 0 ∕E L = 1 + Γsin 2 ( )∕2. As could follow from the crack displacement-stress analysis, E L is independent of the cohesion term. The dependence on y only comes in the expression for the yield stress L y . Typically, for a crack oriented at 45 • to 1 , E 0 /E L ≈ 1 + 1.6Γ. Taking a crack density Γ = 0.5, E L ≈ 0.55E 0 ; that is, the cohesive sliding on the array of cracks causes a 55% modulus deficit relative to the uncracked solid, which is about twice of that observed for the frictional sliding model with the same crack density and a friction coefficient = 0.5.

Unloading
The derivation steps detailed above for loading as the same during unloading, except that for unloading the incremental crack displacement-stress relation (10) must be used instead of (6). The effective Young's modulus E U during unloading of a rock containing a single array of frictional cracks is then where U y is given by Equation 9. In the purely frictional case ( y = 0), typically, taking = 0.5 and = 45 • , E 0 /E U ≈ 1 + 2.4Γ. For a crack density Γ = 0.5, E L ≈ 0.5E 0 ; that is, the frictional backsliding on the array of cracks causes a 50% modulus deficit relative to the uncracked solid, which is greater than the modulus deficit caused by crack sliding during loading, as analyzed in section 2. In the purely cohesive case ( = 0), as could follow from the crack displacement-stress analysis, the effective Young's modulus during unloading is the same as during loading, except that it applies to a different stress range.

Reloading
Consider now reloading of a rock containing an array of frictional cracks that has been previously subjected to a load-unload cycle to a maximum stress ⋆ 1 . The behavior of all cracks during reloading will be identical since they all have the same orientation. Following the detailed analysis of the behavior of a frictional crack during reloading given in section 2.3, and similar derivation steps as those detailed above, the effective Young's modulus during reloading will be (i) If ⋆ 1 < L y (no sliding during the previous load-unload cycle), where L y is given by Equation 5. E RL is then simply given by Equation 20. (ii) If L y < ⋆ 1 < RL y (forward sliding only during the previous load-unload cycle), where RL y is given by Equation 14. (iii) If ⋆ 1 ≥ RL y (backsliding during the previous load-unload cycle),

Body With Randomly Oriented Cracks 3.2.1. Loading
Now imagine that the rock contains a distribution of N frictional cracks per representative area A, each with the same length and with their orientation angles uniformly distributed. Friction on each crack is described by the constitutive law (2). By symmetry, we need only to consider the range of values 0 ≤ ≤ /2. We invoke again the "no-interaction" approximation to add individual crack contributions. However, when cracks are randomly oriented, the behavior of the rock during loading is more complicated than that with aligned cracks, because, as shown in Equation 5, each crack will start sliding at its "own" orientation-dependent yield stress. To calculate the effective Young's modulus, a modified formulation of Equation 20 must therefore be given to account for angular limits of sliding activity. Figure 3 shows the crack yield stress during loading, L y (Equation 5, normalized to 2 ) as function of the crack orientation, in the purely frictional case ( y = 0) for a friction coefficient = 0.5 (Figure 3a) and in the purely cohesive case ( = 0) for a cohesion stress y = 2 (Figure 3b). Analysis of Equation 5 for crack endowed with both frictional and cohesive resistance reveals that L y is a minimum at = (1∕2)[ ∕2 + tan −1 ( )], which is the optimal orientation for crack sliding, and such minimum is equal to which corresponds to the yield stress of the rock with randomly oriented frictional cracks, during loading.
In the purely frictional case ( y = 0), the expressions for the optimal angle for sliding and the yield stress as given above both match those previously derived using the "sliding crack" model (e.g., Ashby & Sammis, 1990;Kachanov, 1982a). In the purely cohesive case ( = 0), the optimal orientation for crack sliding is independent of y and equal to = 45 • .

10.1029/2019JB018970
An equivalent but complementary way of considering the schematic plot of Figure 3 is to identify the yield stress on the vertical axis to a given applied stress to the rock, 1 . At such stress, all individual cracks whose "own" yield stress is below 1 have started sliding, and those whose "own" yield stress is above 1 will require a higher stress to start sliding. Hence, the range of sliding crack orientations is stress dependent and determined by the intersection of a given stress with the yield envelope as shown in Figures 3a and 3b for the purely frictional and cohesive cases, respectively. Recalling previous analysis, crack sliding can occur if the effective sliding stress L eff is positive, where L eff is given by Equation 4. Examination of the range of crack orientations for which L eff > 0 at a given stress reveals that, during loading, sliding can occur between two critical angles ( L 1 , L 2 ) explicitly given by Note that, in the purely frictional case ( y = 0), such results are equivalent to those obtained using a "sliding on a plane of weakness" model (Jaeger et al., 2007). The critical angles for sliding are stress dependent, and the proportion of cracks that can slide increases with increasing stress (Figures 3a and 3b). It can be shown by manipulation of expressions (25) and (26) that tan −1 ( ) ≤ L 1 ≤ L 2 ≤ ∕2 and that, as 1 → ∞, L 1 → tan −1 ( ) and L 2 → ∕2. No sliding occurs in the crack orientation range 0 < < tan −1 ( ) regardless of the applied stress. The value tan −1 ( ) is the "friction angle" (Walsh, 1965). In the purely cohesive case ( = 0), as 1 → ∞, L 1 → 0, and L 2 → ∕2, hence, the entire range of crack orientations can, in principle, slide under sufficiently elevated stress.
Having determined the angular limits for sliding activity and the rock yield stress during loading, if the rock contains randomly oriented frictional cracks, we must integrate the right-hand term of Equation 19 for cracks that are sliding, that is, L 1 ≤ ≤ L 2 . Hence, recalling that we consider the range 0 ≤ ≤ /2, for 1 ≥ L y Equation 19 becomes where we make use of the assumption that all cracks have the same length. After integration of Equation 27, we find that during loading the Young's modulus of a rock containing randomly oriented frictional cracks is where L y is given by Equation 24. That the inverse of the effective Young's modulus depends linearly on the crack density Γ in Equation 28 is intrinsically due to our formulation of the problem in terms of the additional compliance caused by cracks in Equation 18, and that we sum the contributions of sliding cracks assuming no interactions.
The "maximum modulus deficit" caused by crack frictional sliding can be estimated by taking the limit of Equation 28 as 1 → ∞ and compared in the purely frictional and cohesive cases. In the purely frictional case ( y = 0), according to previous analysis, sliding occurs over tan −1 ( ) ≤ ≤ ∕2, and from Equation 28 E 0 ∕E L → 1 + (Γ∕2)[ ∕2 − tan −1 ( ) − ∕(1 + 2 )] as 1 → ∞. Taking a friction coefficient = 0.5, E 0 /E L ≈ 1 + 0.4Γ. For a crack density Γ = 0.5, the modulus deficit relative to the uncracked solid is about 25%. In the purely cohesive case ( = 0), as 1 → ∞ sliding occurs over 0 < < /2, and from Equation 28 E 0 /E L → 1 + Γ/4 ≈ 1 + 0.8Γ. Such comparison between the purely frictional and cohesive cases shows immediately that, as a rule of thumb, in order to produce an equivalent modulus deficit the purely frictional model will require a crack density about twice as that of the purely cohesive model.

Unloading
During unloading, as the load decreases each crack will start backsliding at its "own" orientation-dependent yield stress. Figure 3 shows the crack yield stress during unloading, U y (Equation 9, normalized to 2 ) as function of the crack orientation, considering a rock previously loaded to a maximum stress ⋆ 1 = 10 2 , in the purely frictional case ( y = 0) for a friction coefficient = 0.5 (Figure 3a) and in the purely cohesive case ( = 0) for a cohesion stress y = 2 (Figure 3b). The optimal crack orientation for backsliding is that for which U y is a maximum. Applying the superposition principle as described in section 2.1 to Equation 20b of Lawn and Marshall (1998), such maximum occurs at an angle = tan −1 [ and is equal to the unloading yield stress of a rock containing randomly oriented frictional cracks, Both the optimal orientation for backsliding and the unloading yield stress do not only depend on the friction coefficient and cohesion stress y , but also on the maximum stress experienced during previous loading, ⋆ 1 . In the purely cohesive case ( = 0), the optimal crack orientation for backsliding is, as for loading, independent of the model parameters and equal to = 45 • .
Following a similar analysis to that done for loading above, examination of the criterion U eff > 0 in Equation 8 reveals that, at given stress during unloading, backsliding occurs between two stress-dependent critical angles ( U 1 , U 2 ) given by By analogy with the passing from Equation 19 to Equation 27 for loading, but this time making use of Equation 21, after integration we find that the effective Young's modulus during unloading of a rock containing randomly oriented frictional "backsliding" cracks is given by where U y is given by Equation 29.

Reloading
As in section 3.1.3 consider reloading after previous application of a load-unload cycle to a maximum stress ⋆ 1 , but this time for a rock containing randomly oriented cracks. Based on the analysis given in section 2.3, at the onset of reloading the cracks can be exhaustively classified in three "families": cracks for which no sliding occurred, cracks for which only forward sliding occurred, and cracks for which backsliding occurred-the proportion of which is dictated by the maximum stress during previous loading.
For the crack family for which backsliding occurred, Figure 3a shows the orientation-dependent crack yield stress during reloading, RL y (Equation 14, normalized to 2 ), in the purely frictional case ( y = 0) for a friction coefficient = 0.5 (Figure 3a) and in the purely cohesive case ( = 0) for a cohesion stress y = 2 (Figure 3b). By finding the minimum of RL y , the optimal orientation for sliding upon reloading is, as for loading, = (1∕2)[ ∕2 + tan −1 ( )], and the reloading yield stress of a rock containing randomly oriented cracks is As can be seen by comparing Equations 33 and 24, in a rock containing randomly oriented cracks the yield stress associated with the reactivation of sliding during reloading is always higher than that during a previous 10.1029/2019JB018970 initial loading (as shown in Figures 3a and 3b for the purely frictional and cohesive cases, respectively). Following previous analysis, examination of the criterion RL eff > 0 in Equation 13 shows that, at a given stress, the cracks that are reactivated by sliding during reloading are oriented between two critical angles ( RL 1 , RL 2 ) given by Hence, based on all previous analysis and as illustrated, respectively, in the diagrams shown in Figures 3a  and 3b for the purely frictional and cohesive cases, at the onset of reloading the three possible crack configurations are precisely defined as follows. The cracks that did not slide during the previous cycle are those )) are given by Equations 25 and 26, here taken at the maximum stress ⋆ 1 . The cracks that slid but did not backslide during the previous cycle are then those for which L ( 1 = 2 )) are given by Equations 30 and 31, taken at the end of unloading ( 1 = 2 ). It can be shown, from comparison of expressions (25) and (30), and (26) and (31), respectively, that this family of cracks always exists, as L 1 ( 1 = ⋆ 1 ) < U 1 ( 1 = 2 ) and U 2 ( 1 = 2 ) < L 2 ( 1 = ⋆ 1 ) are both strict inequalities. Finally, cracks for which U 1 ( 1 = 2 ) < < U 2 ( 1 = 2 ) have previously backslid. Continuing on previous analysis, and providing that the previous load-unload cycle was performed at a sufficiently high stress so that backsliding occurred, when the stress is increased the rock has two distinct yield stresses during reloading. The first yield stress, RL y given by (33), is associated with the onset of sliding on cracks that previously backslid. The proportion of these sliding cracks increases with stress and is given by the critical angles for reloading (34) and (35). The second yield stress, ⋆ 1 , corresponds to the onset of sliding on all cracks that slid during the previous load cycle. By comparison of expressions (30) and (34), and (31) and (35), respectively, RL 1 ( 1 = ⋆ 1 ) = U 1 ( 1 = 2 ) and RL 2 ( 1 = ⋆ 1 ) = U 2 ( 1 = 2 ). This implies that, when the stress 1 reaches again ⋆ 1 , all cracks that previously backslid have been reactivated by sliding during reloading. Hence, when 1 = ⋆ 1 , the populations of "sliding only" and "backsliding" cracks merge. Sliding is activated for all cracks for which L 1 ( 1 = ⋆ 1 ) < < L 1 ( 1 = ⋆ 1 ), and, in agreement with the analysis given for the single crack case in section 2.4, the structure of the rock is entirely reset.
For randomly oriented cracks, by analogy with previous derivations for loading and unloading, and making use of Equations 22 and 23 and the analysis given above, the effective Young's modulus during reloading will be (i) If ⋆ 1 < L y (no sliding during the previous load-unload cycle), where L y is given by Equation 24. E RL is then simply given by Equation 28. (ii) If L y < ⋆ 1 < RL y (forward sliding only during the previous load-unload cycle), where RL y is given by Equation 33. (iii) If ⋆ 1 ≥ RL y (backsliding during the previous load-unload cycle),

Effect of the Model Parameters
For a rock containing cracks endowed with both frictional and cohesive resistance and that are either aligned (section 3.1) or randomly oriented (section 3.2), closed-form expressions have been derived for the evolution of the elastic modulus with the applied state of stress, for the loading, unloading, and reloading regimes. As the elastic modulus is stress dependent, the strain can be calculated by integrating the relation d 1 = d 1 /E. For a body with aligned cracks, the effective Young's modulus during loading, unloading, and reloading is only stress dependent through the associated yield stresses demarcating the elastic and inelastic regimes. Both such regimes are linear with stress; hence, analytical expressions are easily obtained for each of the linear segments of the stress-strain curve. In the case of a body with randomly oriented cracks, the effective Young's modulus during loading, unloading, and reloading is only stress dependent through the associated critical angles which account for the proportion of cracks that are sliding. That the proportion of such sliding cracks gradually "stretches" with stress ( Figure 3) causes nonlinearity in the stress-strain curve in the inelastic regime. However, by looking at the structure of the closed-form expressions for the elastic modulus, for a random distribution of crack orientation, it is expected that at sufficiently high stress the inelastic regime must eventually become linear with stress, when all sliding cracks have been activated. Nevertheless, except for the elastic regime, for randomly oriented cracks the strain must be calculated by numerical integration.
For cyclic loading of rock, the modulus and the following stress-strain curve depend only on the following parameters: E 0 , the Young's modulus of the uncracked rock; Γ, the crack density; , the friction coefficient acting on the crack faces; and y , the cohesive yield stress for the cohesive sliding model; , the crack orientation, is the additional parameter for the case of aligned cracks. E 0 simply enters as a "normalizing factor" in all closed-form expressions giving the effective Young's modulus E; hence, the dependence on E 0 is trivial. For aligned cracks, the way the crack orientation enters the model is already captured in the crack displacement-stress analysis given above and shown in Figure 2.
For the purely frictional sliding model ( y = 0), Figure 4a shows the stress-strain curves during two successive load-unload cycles to a maximum stress ⋆ 1 = 10 2 , for two values of the crack density (Γ = 0.5; Γ = 1) at fixed friction coefficient ( = 0.5), for the case of cracks aligned at = 45 • to 1 . Figure 4b shows the equivalent plot for three values of the friction coefficient ( = 0; = 0.3; = 0.6) at fixed crack density (Γ = 1). Figures 4c and 4d show the equivalent plots of, respectively, Figures 4a and 4b for the case of randomly oriented cracks. From our definitions of normalized stress as 1 / 2 and normalized modulus as E/E 0 , it follows that the normalized strain can be defined as (E 0 / 2 ) 1 . Accordingly, a slope of 1 in Figure 4 indicates a normalized modulus of 1, that is, purely elastic behavior (E = E 0 ).
Because we use the noninteraction approximation, the effective compliance E 0 /E is a linear function of the crack density Γ. The resulting effect of increasing the crack density is a more pronounced compliance in the inelastic regime with the yield point invariant (during both loading and unloading segments), and a larger hysteresis loop with more "permanent" inelastic strain at the end of unloading (Figures 4a and 4c). As can be physically expected, an increasing friction coefficient results in a less pronounced compliance in the inelastic regime for both loading and unloading, with a higher (resp. lower) yield point during loading (resp. unloading) leading to a greater proportion of elastic behavior relative to inelastic behavior and, generally, a less pronounced hysteresis loop (Figures 4b and 4d). In the case of aligned cracks, using the analysis given in section 2.3 and Equation 16, it can be shown that the permanent strain at the end of unloading is exactly equal to Γ sin(2 ) and thus increases linearly with the friction coefficient, as shown in Figure 4b. In the case of randomly oriented cracks, however, the way the residual inelastic strain varies with the friction coefficient is less trivial (Figure 4d). The case = 0 corresponds to perfectly reversible linear behavior, with a pronounced modulus deficit relative to the uncracked solid (Figures 4b and 4d). In the case of randomly oriented cracks, by manipulation of Equation 28 combined with (25) and (26) in the limit → 0, E 0 /E L → 1 + ( /4)Γ.
The equivalent plot of Figure 4 for the purely cohesive sliding model ( = 0) is shown in Figure 5. The fixed values of the cohesive yield stress when the crack density is varied are y ∕ 2 = 1 (Figures 5a and 5c), and the three representative values of y at fixed crack density are y ∕ 2 = 0, y ∕ 2 = 1, and y ∕ 2 = 2 (Figures 5b and 5d). The resulting effect of increasing the crack density (Figures 5a and 5c) is the same as that described above for the purely frictional model. Increasing the cohesive yield stress results in a higher (resp. lower) yield point during loading (resp. unloading) but does not affect the compliance in the inelastic regime (Figures 5b and 5d). For both aligned and randomly oriented cracks, the permanent strain at the end of unloading increases with the cohesive yield stress. In the case of aligned cracks, using the analysis given in section 2.3 and Equation 16, the permanent strain at the end of unloading is found to be equal to Γ y sin(2 ). The case y = 0 corresponds to perfectly reversible linear behavior, with a pronounced modulus deficit relative to the uncracked solid (Figures 5b and 5d). In the case of randomly oriented cracks, by manipulation of Equation 28 combined with (25) and (26) in the limit y → 0, E 0 /E L → 1 + ( /4)Γ. As could be physically expected, such limit is the same than that obtained for the equivalent purely frictional sliding model as → 0.
As previously analyzed, for both purely frictional and cohesive sliding models, respectively, Figures 4 and 5 illustrate that, when the maximum stress ⋆ 1 is reached during reloading the rock strain is the same as that at ⋆ 1 during previous loading; hence, the stress-strain behavior is reversible upon reloading. Including a random distribution of cracks results in a "smoother" transition between the elastic and inelastic regimes in the stress-strain behavior (around the yield point), which corresponds to a progressive, stress-dependent activation of sliding or backsliding on favorably oriented cracks. The observed difference in the stress-strain path between the cases of aligned and randomly oriented cracks is also due to the fact that the contribution of the crack sliding displacement to the axial strain is orientation dependent, as shown by Equation 16. Finally, note that the more pronounced compliance and hysteresis for the case of aligned cracks directly result from the fact that the selected crack orientation ( = 45 • ) for sliding is near optimal in the frictional case, and optimal in the cohesive case.

Comparison of the Stress-Strain Behavior for the Purely Frictional and Cohesive Sliding Models
A comparison of the stress-strain behavior for the purely frictional and cohesive models is shown in Figure 6, taking the same fixed value of crack density Γ = 1. A fixed value of the coefficient of friction, = 0.5, and its "equivalent" cohesive yield stress, y ∕ 2 = 0.5 (n.b. "equivalent" in the way that, in the case of Figure 5. Stress-strain curves in the purely cohesive case ( = 0) for two successive load-unload cycles at a maximum stress ⋆ 1 = 10 2 , for aligned cracks oriented at = 45 • to 1 (a, b), and for randomly oriented cracks (c, d). Showing the effect of crack density (Γ = 0.5; Γ = 1) at fixed cohesion ( y ∕ 2 = 1) (a, c), and the effect of cohesion ( y ∕ 2 = 0; y ∕ 2 = 1; y ∕ 2 = 2) at fixed crack density (Γ = 1) (b, d). Solid lines: loading segments; dashed lines: unloading segments. aligned cracks, as demonstrated above the same permanent strain is produced for both models), is used for the purely frictional and cohesive models, respectively. The additional parameter for the case of aligned cracks (a, c) is the crack orientation taken as = 45 • . As the detailed features of the stress-strain behavior have already been described above, and notably the effect of including a crack distribution, here we highlight the main differences demarcating the purely frictional and cohesive models. (i) At a given crack density, in the inelastic regime the purely cohesive model predicts a greater modulus deficit relative to the uncracked solid than the purely frictional sliding model. The inelastic compliance depends on the friction coefficient for the purely frictional model but remains independent of y for the purely cohesive model. (ii) The purely cohesive model predicts the same inelastic compliance during both loading and unloading, whereas the purely frictional model predicts a greater inelastic compliance during unloading than during loading.

Cyclic Loading of a Rock Containing Multiple Cracks: Dissipated Strain Energy
The dissipated strain energy, W, is defined as the difference between the amount of energy that is recovered on unloading and the amount of energy input on subsequent reloading; that is, W represents the area between the reloading and unloading curve at a given stress. Defined by analogy with internal friction for materials exhibiting "static hysteresis" (Nowick, 1954), the quantity W is thus independent of the number of cycles. W is calculated per unit volume by integration of the stress-strain relations, as W = ∫ 1 d 1 , along the unloading-reloading contour. Considering the features of the crack displacement-stress and stress-strain behavior in the present crack sliding model, in order to dissipate energy upon unload-reload backsliding must necessarily occur. According to previous analysis, this condition is expressed as ⋆ 1 > RL y , where Figure 6. Stress-strain curves for two successive load-unload cycles at a maximum stress ⋆ 1 = 10 2 , for the frictional (a, b) and cohesive (c, d) crack sliding model, taking a crack density Γ = 1 for both models. (a, c) Aligned cracks oriented at = 45 • to 1 ; (b, d) randomly oriented cracks. Taking a friction coefficient = 0.5 for the frictional model, and an "equivalent" cohesive yield stress y ∕ 2 = 0.5 for the cohesive model. Solid lines: loading segments; dashed lines: unloading segments. Squares: yield stresses.

Purely Frictional Sliding
In the case of purely frictional sliding ( y = 0), for a rock containing aligned cracks, analytical expressions can be obtained for the dissipated energy upon unload-reload: where the coefficients (M L , M U ) are expressed as in which it is immediately seen that W is quadratic in ⋆ 1 . For randomly oriented cracks, as for the stress-strain curve the dissipated energy W must be calculated by numerical integration. A comparison of the dissipated energy versus stress behavior for the purely frictional and cohesive models is shown in Figure 7.

Purely Cohesive Sliding
In the case of purely cohesive sliding ( = 0), for a rock containing aligned cracks, analytical expressions can also be obtained for the dissipated energy upon unload-reload: in which it is immediately seen that W is linear in ⋆ 1 . Note that this expression is very comparable to that obtained for a proportional loading (Hansen et al., 2020). . Dissipated strain energy upon unload-reload as function of maximum stress per cycle, using a crack density Γ = 1 for both the purely frictional (black curves) and purely cohesive (gray curves) crack sliding models, and a friction coefficient = 0.5 and cohesive yield stress y ∕ 2 = 0.5, respectively. Full lines: randomly oriented cracks; dashed lines: array of aligned crack at = 45 • to 1 .
For randomly oriented cracks, as for the stress-strain curve the dissipated energy W must be calculated by numerical integration.

Comparison of Purely Frictional and Cohesive Sliding Models
A comparison of the dissipated strain energy versus stress behavior for the purely frictional and cohesive models is shown in Figure 7, taking for each model the same parameter values as in the comparison illustrated in Figure 6 for the stress-strain behavior. For the purely frictional model ( y = 0), the dissipated strain energy W increases with stress and, for values of stress typically greater than 10 times the confining pressure, W ∝ ( ⋆ 1 ) 2 . For the purely cohesive sliding model ( = 0), W also increases with stress but, contrary to the frictional model, in the high stress limit W ∝ ⋆ 1 . Such quadratic and linear dependence of W upon ⋆ 1 for the frictional and cohesive models, which is obvious from the form taken by Equation 37 and 40 for aligned frictional cracks, respectively, is also observed for randomly oriented cracks regardless of the values taken by the model parameters.
For both purely frictional and cohesive models, including a random distribution of crack orientation essentially results in offset in the dependence of W with stress.
Hence, in addition to the differences highlighted above between the purely frictional and cohesive models in terms of stress-strain behavior, in terms of dissipated strain energy we conclude that the striking difference is the quadratic and linear dependence of W upon ⋆ 1 for the purely frictional and cohesive models, respectively, at moderate to high stress. At low stress, W is greater for the purely cohesive model than for the purely frictional model; however, this trend reverses at moderate to high stress due to the high-order dependence of W upon ⋆ 1 for the purely frictional model.

Application of Model to Experimental Data
We now test our sliding crack model on stress-strain data obtained during cyclic loading triaxial experiments on isotropic polycrystalline antigorite specimens. The choice of such rock type and of the applied range of confining pressure and compressive stress is intended to satisfy the model hypothesis (a-d, see section 1; also, see section 6). In addition to fitting the model to the cyclic loading data obtained by David et al. (2018) on an isotropic polycrystalline antigorite sample at room temperature and 150 MPa confining pressure, we have also carried out a set of cyclic loading, triaxial compression experiments in the Griggs apparatus at the Department of Earth, Environmental and Planetary Sciences at Brown University (Providence). The rock selected is the same material as the "isotropic" block of "Vermont antigorite serpentinite" described in detail by David et al. (2018). The rock is fine grained, nearly pure antigorite (95%) with minor amount of magnetite and magnesite, and is elastically isotropic.
Three cylindrical core specimens were precision ground to 12.7 mm length and 6.17 mm diameter. One thin Ni disk (thickness of 0.2 mm) was placed at each end between the rock sample and alumina (Al 2 O 3 ) pistons. This assembly was enclosed within a thin-walled (0.254 mm wall thickness) silver jacket. For the two experiments conducted at 400 • C and 500 • C, a solid salt (NaCl) assembly was used as a confining medium both inside and outside the graphite furnace, and sample temperature was monitored using a Pt-Pt10%Rh thermocouple (see Figure 1b of Holyoke & Kronenberg, 2010, for additional details on the sample assembly). For the room temperature experiment, an all-lead assembly was used as the confining medium. The specimen assembly was then loaded in a Griggs apparatus (see, e.g., Holyoke & Kronenberg, 2010) for cyclic loading axial compression deformation experiments at room temperature, 400 • C and 500 • C, and at a confining pressure of 1 GPa (accuracy of 0.05 GPa). The 400 • C and 500 • C samples were brought to pressure and temperature over 3-4 hr, to 100 • C at 250 MPa, to 200 • C at 400 MPa, to 300 • C at 600 MPa, and then to the desired pressure and temperature of the experiment. Following a "prehit" stage, samples were cyclically loaded at a constant strain rate of 10 −5 s −1 , gradually increasing the maximum stress per cycle to about 90-95% of 10.1029/2019JB018970 sample strength, the latter which was determined from additional axial deformation experiments at larger strains (not reported here). At the end of the experiment, pressure and temperature were reduced slowly.
The axial shortening was measured with an external linear variable differential transformer (LVDT) and corrected from machine compliance to measure rock specimen axial strain. The machine compliance was refined by an additional calibration experiment performed on a WC sample (in a solid salt assembly) that was axially loaded in its elastic range at 400 • C and 500 • C, and found to be equal to 9.0 μm·kN −1 . The axial load was measured externally by a load cell that was recalibrated for this series of experiments, converted to axial stress by dividing by the sample diameter, and subtracted by the confining pressure. The obtained "raw differential stress" was then corrected for the apparatus frictional sliding forces resulting from advancement of the axial 1 piston. The hit-point was found following the method described in Holyoke and Kronenberg (2010) (see Figure 3a of that publication) as the intersection between the run-in through lead (initially present between the alumina top piston and the WC loading piston) and the elastic loading of sample, from which friction was directly estimated. When piston direction is reversed at the onset of unloading, a portion of stress-strain behavior was removed corresponding to the need for the piston to overcome twice the frictional sliding force, as well as a small displacement offset associated with rotation of LVDT bracked during the change in direction of the motor used to generate the load. Although friction has been reported to increase with axial displacement (Holyoke & Kronenberg, 2010), a thorough examination of the cyclic loading data at the onset of unloading indicated that no such correction were needed due to the overall small displacements involved. The differential stress was then finally corrected for changes in sample area, by taking a Poisson's ratio equal to 0.28 from David et al. (2018). The obtained differential stress on sample is considered to be known to ±50 MPa.
The strategy and assumptions in fitting the data were as follows. For both the purely frictional ( y = 0) and purely cohesive ( = 0) models, since the rock is elastically isotropic, a random orientation of cracks was assumed. A pressure-independent value of the Young's modulus of the uncracked solid, E 0 , was first determined at each temperature (20 • C, 400 • C, and 500 • C) by fitting the linear elastic portion of the stress-strain curves at low stress (preyield). Since the rock is isotropic, for both purely frictional and cohesive models, a single value of the crack density Γ was used for each model at all temperatures and pressures. Γ is an intrinsic property of the rock that should, to first order, not depend on pressure, under the assumption that pressure-induced crack closure is reasonably complete in antigorite at 150 MPa (David et al., 2018), and under the model assumption that no crack propagation occurs. It is also assumed that Γ does not vary with temperature, as the process of thermal cracking and generally the presence of open microcracks must be inhibited by the substantial confining pressures applied during the experiments. For the purely frictional and cohesive models, respectively, it seems reasonable to assume that both the friction coefficient and cohesive yield stress y only depend on temperature, since temperature is likely to affect the nature of physical bonds between the two crack faces. Note that some portions of the stress-strain curves were discarded from the fit based on the following reasons. The first situation is an excessively pronounced downward inflexion of the stress-strain curve during the loading cycle at the highest stress, which is possibly caused by the onset of crack propagation and thus not captured by the model. The second reason is that a substantial portion of the unloading segments (particularly at the onset of unloading) for the data obtained in the Griggs' apparatus yield unphysically large Young's moduli that are well above that extracted from the elastic portion of the loading segments; modification of the data processing strategy to extract more physically plausible stress-strain data during unloading appears to be an excessively subjective and ad hoc process. Hence, for those stress-strain data discarded from the fit, the model curves are predictions.
By least squares inversion, the best fit of the data for each model is obtained using the fitting parameters given in Table 2 for the entire data set and shown for each individual cycle in Figure 8. The posterior (unnormalized) probability density function in the (Γ, ) and the (Γ, y ) spaces for the purely frictional and cohesive models, respectively, is shown in Figure 9, at each temperature. Considering the highly nonlinear and hysteretic shape of the stress-strain curve, it can be claimed that both models do an overall good job in fitting the multiple loading-reloading curves at all conditions of pressure and temperature. At room temperature, at 150 MPa confining pressure the data are better fit by the purely frictional model than by the purely cohesive model, and the purely frictional model does a good job in capturing the unloading curves, whereas at 1 GPa confining pressure the data are better fit by the purely cohesive model than by the purely frictional model, and the purely cohesive model provides a good fit to the unloading curves for the first two cycles. At 400 • C and 500 • C, both models provide very good and essentially identical fit to all loading-reloading curves except the high stress portion of the last two cycles at 500 • C. Hence, overall it can be said that the purely frictional model seems to do a better job in capturing the data at lower-pressure (150 MPa) data, whereas the purely cohesive model is only able to describe the data at higher pressure (1 GPa).
The model-independent Young's modulus of the uncracked rock decreases with temperature, from 103 GPa at room temperature to 83 GPa at 500 • C. As might have been expected from previous analysis, the crack density obtained from the inversion is greater for the purely frictional model (0.93) than for the purely cohesive model (0.65). The inverted friction coefficient decreases with temperature, from 0.41 at room temperature to 0.16 at 500 • C. The temperature dependence of the cohesive yield stress is less clear than that of the friction coefficient, but y broadly decreases from 243 MPa at room temperature to 182 MPa at 500 • C. For the purely frictional and cohesive models, respectively, the inverted values of the friction coefficient and cohesive yield stress are much better constrained than that of the crack density ( Figure 9).

Physical Significance of the Inverted Model Parameters for Antigorite
The values of the model parameters used in fitting the cyclic loading stress-strain data on antigorite (Figure 8) all seem physically reasonable. The temperature-dependent values of E 0 , the Young's modulus of the uncracked rock, are given in Table 2. At room temperature, Bezacier et al. (2010) report E 0 = 97 GPa at room pressure using aggregate averages from single-crystal elasticity data, and David et al. (2018) report E 0 = 93 GPa at 150 MPa confining pressure (i.e., for the data set shown in Figures 8a-8d) using ultrasonic wave velocity measurements. Both values of E 0 are consistent with the value found here, E 0 = 103 GPa. David et al. (2019) recently reported a temperature dependence of shear modulus of about −17 MPa·K −1 using low-frequency torsional oscillation tests, on exactly the same material as that used here for the cyclic loading data. This gives a decrease in shear modulus of about 8.2 GPa in the 20-500 • C range. If a temperature-independent Poisson's ratio of 0.28 is assumed (David et al., 2019), using elastic constant relationships the corresponding decrease in Young's modulus over the 20-500 • C is 21 GPa, which is consistent with what we found here. The values for the friction coefficient between the crack faces (Table 2) can be compared to representative values obtained from deformation experiments on intact or precut antigorite specimens, or antigorite gouge friction experiments. At room temperature, reported values of vary between, for example, = 0.77 (Dengo & Logan, 1981), = 0.5-0.85 (Reinen et al., 1994), and = 0.34 (Escartín et al., 1997). In the 25-200 • C range, Moore et al. (1997) reported = 0.4-0.6. More recently, Chernak and Hirth (2010) reported = 0.15 at 550 • C, and Proctor and Hirth (2016) have measured = 0.23, = 0.13, and = 0.07 at 300 • C, 400 • C, and 500 • C, respectively. Hence, the range of temperature-dependent values of that we obtain using the purely frictional model is consistent with previously published experimental data.
Comparison of the obtained values for the cohesive yield stress with experimental data is difficult due to the scarcity of the available measurements, but also because the physical origin of the cohesion term could arise from several mechanisms (Lawn & Marshall, 1998). Nevertheless, we note that the cohesion y is about 0.7% of the antigorite's shear modulus at room temperature. Hence, the ratio is in the same order of magnitude as the ratio of the "lattice friction" or "Peierls stress" to shear modulus, assuming that such a ratio in antigorite is comparable to other existing measurements on silicates or carbonates (Nabarro, 1997) or computational estimates (Skelton & Walker, 2019).
The best fitting crack densities are 0.93 and 0.65 for the purely frictional and cohesive models, respectively. Such values are reasonable and do not fall far outside the limit of validity of the no-interaction approximation (Guéguen & Schubnel, 2003;Kachanov, 2007). However, they are to some extent artificially high having Figure 8. Fits of the entire stress-strain data set at all temperatures, using a single crack density, and temperature-dependent values of and y for the purely frictional (red) and purely cohesive (blue) crack sliding models, respectively, and the same temperature-dependent Young's modulus E 0 (T) for the two models. Fitting parameters are given in Table 2. Each cycle is shown individually, by increasing maximum stress. (a-d) P c = 150 MPa and room temperature (data from David et al., 2018); (e-g) P c = 1,000 MPa and room temperature (data from this study); (h-m) P c = 1,000 MPa and T = 400 • C (data from this study); (n-s) P c = 1 000 MPa and T = 500 • C (data from this study). Experimental data used (resp. discarded) for fit: black (resp. gray) full circles. Model curves: load (full line); unload (dashed line). Gray dashed line: purely elastic stress-strain curve. neglected stress-field interactions between nearby cracks. If we had instead used the differential effective medium approximation, the rescaled crack density would be (1∕ ) ln(1 + Γ) (David et al., 2012), for example, Γ = 0.44 and 0.35 for the purely frictional and cohesive models, respectively. In general, it is expected that using an effective medium theory for taking crack interactions into account would mostly result in rescaling the crack density but would not fundamentally alter the nature and the shape of the evolution of the physical properties with stress, here the Young's modulus. The contentious debate of which effective medium theory best accounts for interactions between sliding cracks is beyond the scope of our work.
Finally, we note that the fitted parameters for the purely frictional sliding model at room temperature are in good agreement with those inverted from indentation data on antigorite single crystals (Hansen et al., 2020) from exactly the same rock material as that used in the experiments of Figure 8. Hansen et al. (2020) applied both purely frictional and cohesive sliding models to the case of a single array of cracks, but under conditions of proportional loading that are more representative of self-confined indentation experiments. They found an uncracked Young's modulus of 97 GPa; for the purely frictional sliding model, a friction coefficient of 0.5 and crack densities in the range 0.2-0.6 depending on crystals; and for the purely cohesive model, a cohesion of 150 MPa and crack densities in the range 0.1-0.4 depending on crystals.

Transition in Sliding Mechanism With Increasing Pressure and Temperature
A striking outcome of the application of the crack sliding model to the data is that the stress-strain behavior at all pressure and temperatures in antigorite can be reasonably described by a coefficient of friction, , that decreases with temperature ( Table 2). As noted above, the temperature-dependent values of inverted by utilizing the purely frictional model are in good agreement with previous estimates from friction experiments at high strain. That decreases with temperature carries important implications for the rheology of subduction zones. In particular, Chernak and Hirth (2010) and Proctor and Hirth (2016) have speculated that such a decrease in could directly explain the inverted "brittle to ductile" transition, that is, the transition back from distributed to localized deformation behavior, that has been experimentally observed in antigorite above 300 • C (at pressures greater than a few hundreds of MPa).
Another important outcome of the application of the model to the data is that, globally, the low-pressure (150 MPa) data are better captured by the purely frictional model, while the high-pressure (1 GPa) behavior is overall better described by the purely cohesive model at all temperatures. A transition from a frictionally to a cohesively controlled regime of crack sliding with increasing pressure is consistent with the mechanistic interpretation of the brittle-to-plastic transition in the crust (Escartín et al., 1997;Hirth & Tullis, 1994), providing that the nature of cohesive sliding can be attributed to some form of plasticity, such as overcoming a "lattice friction" or "Peierls stress" as briefly discussed above. Of course, additional data at intermediate pressures and temperatures would be required to reinforce such interpretations.

Limitations of the Model and Onset of Crack Propagation
The inability of both purely frictional and cohesive models to fit the loading curve during the stress-strain cycle at the highest stress in some cases, for example, at 1 GPa and room temperature (Figure 8g) or 500 • C (Figure 8s), is likely caused by the onset of crack propagation in the antigorite samples, resulting in additional inelastic compliance that cannot be captured by the model. In antigorite such process occurs in the stress range typically between 90% and 100% of the rock strength. Previous experimental observations have demonstrated that the nature of such crack propagation in antigorite is peculiar compared to that of other rocks. Up to failure, the brittle deformation of antigorite serpentinites is nondilatant at room temperature at least up to 300 MPa confining pressure (Escartín et al., 1997) and not associated with the opening of Mode I or "wing" microcracks (David et al., 2018). Hence, as proposed in these publications, the cracks that propagate in antigorite are Mode II or in-plane "shear" microcracks. Incorporating Mode II propagation in the crack sliding model at high stress is scope for future work.
That "Mode II" is the favored mode for crack propagation in antigorite must reflect that the onset of dilatancy or "Mode I" crack opening would require a substantially large stress. To gain knowledge on what such stress could be, we recall that for a rock containing randomly oriented microcracks, the stress at which "wing cracks" initiate is given by (Nemat-Nasser & Horii, 1982) where K Ic is the Mode I fracture toughness of the material. As previously noted, the first term on the right-hand side of Equation 41 corresponds to our expression (24) for the yield stress of a rock containing randomly oriented cracks, in the purely frictional case ( y = 0). The second term on the right-hand side of Equation 41 gives the additional stress required to initiate Mode I opening past the yield point, which we refer here as the "dilatancy term".
To further test the hypothesis that in antigorite such term may be quantitatively significant, we have conducted direct measurements of K Ic in polycrystalline antigorite using the semicircular bend (SCB) methodology (Kuruppu et al., 2014) but also measurements of the tensile strength, t , by the Brazil disk test (International Society or Rock Mechanics, 1978), both at room temperature. For both the SCB and Brazil tests, sample geometry and experimental setup are as described in Inskip et al. (2018), to the exceptions that the constant displacement rates are 0.3 and 0.5 mm/min for the SCB and Brazil tests, respectively; the ratio of the notch depth to the sample radius is 0.4 (see Figure 6 of that publication) for the SCB test; and the sample diameter and thickness are 40 by 20 mm for the Brazil test. The SCB and Brazil tests were carried out, respectively, on two and three isotropic antigorite polycrystalline samples of exactly the same material as in the experiments of David et al. (2018) and the Griggs experiments presented above. The fracture toughness of antigorite is found to be K Ic = 2.2 MPa·m 1/2 . Such value appears to be in the high range of K Ic measured on crystalline rocks (Zhang, 2002) and slightly greater than those measured on granites, for example, Westerly granite (K Ic = 1.8 MPa·m 1/2 , Meredith & Atkinson, 1985). The experimentally determined value of the tensile strength from Brazil tests is t = 34 MPa. The valuable information gained from this measurement is that the ratio of the tensile strength to the fracture toughness for antigorite is about 15·m −1/2 , which is about 2.2 times the ratio typically observed for most rocks ( t = 6.9K Ic , Zhang, 2002). We recall that, according to Griffith analysis, t = (C∕ √ 2c)K Ic , where C is a dimensionless geometric constant and c is the characteristic crack radius. If we now invoke the assumption that the geometrical constant C does not drastically change between materials, that the ratio t /K Ic in antigorite is 2.2 greater than most rocks implies a characteristic crack size which is 2.2 2 ≈ 5 times smaller than the mean for other rock materials-at least in tension.
Looking at the functional form of the dilatancy term in Equation 41, it is therefore entirely reasonable to expect that, under a relatively high K Ic and a small crack radius c, the dilatancy term is very large for antigorite. We have accordingly reported in Table 3 the calculated contributions of the "yield term" and the "dilatancy term" to the dilatancy stress 1 in Equation 41, as functions of pressure and temperature, using the obtained model values for the friction coefficient at each temperature (Table 2), the experimentally measured K Ic (assumed to be temperature independent), and a characteristic crack length assumed to be equal to the average observed grain size in the antigorite samples, that is, 2c = 10 μm. At all pressures and temperatures, the dilatancy term is greater than 1 GPa, and the total stress required for the onset of dilatancy is larger than the maximum stress in the cyclic loading experiments (Figure 8), which validate our model assumption and data interpretation that, in antigorite, Mode I crack propagation is not a favored mechanism for producing inelastic strain.

Conclusions
We have extended the microphysical model of David et al. (2012) for the effect of sliding cracks on the uniaxial loading-unloading of a rock to a triaxial state of stress and to the entire reloading process. In addition, analysis has been given for a combined "cohesive plus frictional" constitutive crack sliding behavior. Results have been derived for cracks that are either randomly oriented or aligned at a given angle to the maximum compressive stress. For loading, unloading, and reloading, closed-form expressions have been derived for the evolution of the Young's modulus with stress, for the yield stress demarcating the transition between the elastic and inelastic regimes, and for the critical angles for sliding activity in the case of randomly oriented cracks. The dissipated strain energy upon unloading-reloading has also been calculated. The two constitutive "end cases" of purely frictional and cohesive crack sliding are treated separately and compared to experimental data. Both purely frictional and cohesive models provide a good fit to a cyclic loading data set on the triaxial deformation of polycrystalline antigorite under various pressure and temperature conditions: One portion of the data set, at 150 MPa confining pressure and room temperature, was taken from the literature; the rest of the data set at 1 GPa confining pressure and room temperature, 400 • C and 500 • C, was newly acquired in a Griggs deformation apparatus. We found that the stress-strain curves can be inverted from the models assuming that the specimens have the same crack density, temperature-dependent uncracked Young's modulus, and coefficient of friction (for the purely frictional model) or cohesive yield stress (for the purely cohesive model). The values of the parameters are physically reasonable. The stress-strain behavior at increasing temperature can be simply explained by a decreasing friction coefficient. The purely frictional model seems to better capture the low-pressure behavior, while the purely cohesive model could be more adequate at high pressures, at least in the range of pressure of this study.
The ability of the model to characterize actual numerical values of the parameters is limited by a number of assumptions. First, we assume a plate-like specimen under plane-stress conditions. The choice of a two-dimensional analysis was motivated by the greater simplicity of mathematical treatment, but also by the intended scope of incorporating in-plane crack propagation in the model at a later stage, for which calculations are much more approachable in the two-dimensional case. Extending the present model to a three-dimensional analysis would pose no major conceptual and computational problems but would essentially require additional integrations for crack orientations over the azimuthal range as in Kachanov (1982a). Nevertheless, it is expected that results in an axisymmetric three-dimensional case would essentially differ from those in the two-dimensional case by a geometrical factor that should include the Poisson's ratio of the material but should not affect the qualitative evolution of the elastic modulus with stress and the resulting features of the stress-strain curve. We have used the no-interaction approximation for simplicity of the analysis. This rigorous but simplified assumption unavoidably leads to underestimating the effect of sliding cracks on inelastic deformation (Horii & Nemat-Nasser, 1983), and the obtained crack densities should be considered as indicative. However, contrary to purely elastic processes such as wave propagation in a cracked medium, the applicability of effective medium theories to inelastic and dissipative process such as crack sliding (Horii & Nemat-Nasser, 1983) remains unclear. In addition, numerical simulations indicate that, in taking into account the effect of cracks, the no-interaction approximation appears to be accurate up to high crack densities (Kachanov, 2007).
The model can be applied to triaxial deformation experiments (and particularly cyclic loading experiments) on rocks or materials that are initially isotropic, in which one compressive stress is increased above a constant confining stress. The value of the confining pressure should be sufficient so that microcracks are closed and their faces are into contact. The completion of crack closure can be checked by looking at the evolution of elastic wave velocities with increasing pressure, if available. Ongoing crack closure would also result in a "concave-upward" inflexion of the stress-strain curve at the beginning of loading, like that demonstrated by David et al. (2012) during uniaxial loading. The main restriction of the model is that inelastic behavior is purely caused by sliding on cracks; hence, it should be applied to stress ranges or rock types exhibiting no dilatancy or Mode I opening; the latter can be checked from volumetric strain data or elastic wave velocities, for instance. In observing the stress-strain data, the inelastic portions during loading and unloading should be carefully compared as they can distinguish between a frictional-type behavior (greater slope during unloading than during loading) and a cohesive-type behavior (same slope during loading and unloading). In inverting the cyclic loading data here we have made no direct use of the dissipated strain energy, since it is much more precise to invert an entire stress-strain curve than a single value of the dissipated energy for a given stress. Nevertheless, quantitative use of the dissipated strain energy can be convenient in situations where extraction of the stress-strain curve from raw load-displacement data is difficult, as, for instance, in the cyclic loading indentation tests of Hansen et al. (2020). If measured at increasing stress during "internal friction" stress cycles under static or very low-frequency conditions, the dissipated strain energy can also be useful to distinguish between a purely frictional behavior (quadratic in stress) and purely cohesive behavior (linear in stress).

Data Availability Statement
Experimental data are available from the U.K. National Geosciences Data Centre (https://www.bgs.ac.uk/ services/ngdc/) or upon request to the corresponding author.