# A Brittle Constitutive Law for Long-Term Tectonic Modeling Based on Sub-Critical Crack Growth

Léo Petit, Jean-arthur Olive, Alexandre Schubnel, Laetitia Le Pourhiet, Harsha S Bhat

► To cite this version:

Léo Petit, Jean-arthur Olive, Alexandre Schubnel, Laetitia Le Pourhiet, Harsha S Bhat. A Brittle Constitutive Law for Long-Term Tectonic Modeling Based on Sub-Critical Crack Growth. Geochemistry, Geophysics, Geosystems, 2024, 25, 10.1029/2023gc011229. [https://hal.science/hal-04615630](https://hal.science/hal-04615630)

HAL Id: hal-04615630

https://hal.science/hal-04615630

Submitted on 18 Jun 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A Brittle Constitutive Law for Long-Term Tectonic Modeling Based on Sub-Critical Crack Growth

Léo Petit, Jean-arthur Olive, Alexandre Schubnel, Laetitia Le Pourhiet, Harsha S Bhat

► To cite this version:

Léo Petit, Jean-arthur Olive, Alexandre Schubnel, Laetitia Le Pourhiet, Harsha S Bhat. A Brittle Constitutive Law for Long-Term Tectonic Modeling Based on Sub-Critical Crack Growth. Geochemistry, Geophysics, Geosystems, 2024, 25, pp.e2023GC011229. 10.1029/2023gc011229. hal-04615630

HAL Id: hal-04615630

https://hal.science/hal-04615630v1

Submitted on 18 Jun 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Geochemistry, Geophysics, Geosystems®

## RESEARCH ARTICLE

10.1029/2023GC011229

**Special Section:**

Frontiers in lithospheric dynamics: bridging scales through observations, experiments, and computations

**Key Points:**

* New brittle constitutive law describes the onset of faulting in tectonic simulations

* Model is based on sub-critical growth and interaction of micro-cracks

* Laboratory-derived model parameters can be used to model crustal-scale faulting

**Correspondence to:**

J.-A. Olive,
olive@geologie.ens.fr

**Citation:**

Petit, L., Olive, J.-A., Schubnel, A., Le Pourhiet, L., & Bhat, H. S. (2024). A brittle constitutive law for long-term tectonic modeling based on sub-critical crack growth. *Geochemistry, Geophysics, Geosystems, 25*, e2023GC011229. https://doi.org/10.1029/2023GC011229

Received 7 SEP 2023
Accepted 9 FEB 2024

**Author Contributions:**

**Conceptualization:** Léo Petit, Jean-Arthur Olive, Alexandre Schubnel, Laetitia Le Pourhiet

**Data curation:** Léo Petit

**Formal analysis:** Léo Petit, Jean-Arthur Olive, Alexandre Schubnel, Harsha S. Bhat

**Funding acquisition:** Jean-Arthur Olive, Harsha S. Bhat

**Investigation:** Léo Petit, Jean-Arthur Olive, Alexandre Schubnel, Laetitia Le Pourhiet, Harsha S. Bhat

© 2024 The Authors. Geochemistry, Geophysics, Geosystems published by Wiley Periodicals LLC on behalf of American Geophysical Union. This is an open access article under the terms of the Creative Commons Attribution-NonCommercial-NoDerivs License, which permits use and distribution in any medium, provided the original work is properly cited, the use is non-commercial and no modifications or adaptations are made.

## A Brittle Constitutive Law for Long-Term Tectonic Modeling Based on Sub-Critical Crack Growth

<figure><img src="image_2.png" /></figure>

Léo Petit¹ ID, Jean-Arthur Olive¹ ID, Alexandre Schubnel¹ ID, Laetitia Le Pourhiet² ID, and Harsha S. Bhat¹ ID

¹Laboratoire de Géologie, CNRS - École normale supérieure - PSL University, Paris, France, ²Sorbonne Université, ISTEP, Paris, France

**Abstract** Adequate representations of brittle deformation (fracturing and faulting) are essential ingredients of long-term tectonic simulations. Such models commonly rely on Mohr-Coulomb plasticity coupled with prescribed softening of cohesion and/or friction with accumulated plastic strain. This approach captures fundamental properties of brittle failure, but is overly sensitive to empirical softening parameters that cannot be determined experimentally. Here we design a brittle constitutive law that captures key processes of brittle deformation, and can be straightforwardly implemented in standard geodynamic models. In our Sub-Critically-Altered Maxwell (SCAM) flow law, brittle failure begins with the accumulation of distributed brittle damage, which represents the sub-critical lengthening of tensile micro-cracks prompted by slip on pre-existing shear defects. Damage progressively and permanently weakens the rock's elastic moduli, until cracks catastrophically interact and coalesce up to macroscopic failure. The model's micromechanical parameters can be fully calibrated against rock deformation experiments, alleviating the need for ad-hoc softening parameters. Upon implementing the SCAM flow law in 2-D plane strain simulations of rock deformation experiments, we find that it can produce Coulomb-oriented shear bands which originate as damage bands. SCAM models can also be used to extrapolate rock strength from laboratory to tectonic strain rates, and nuance the use of Byerlee's law as an upper bound on lithosphere stresses. We further show that SCAM models can be upscaled to simulate tectonic deformation of a 10-km thick brittle plate over millions of years. These features make the SCAM rheology a promising tool to further investigate the complexity of brittle behavior across scales.

## 1. Introduction

Tectonic plates tend to be almost rigid and primarily deform within narrow boundary zones. In the upper crust (above ~15 km depth), deformation occurs in the brittle regime through nucleation and growth of fractures and faults, which profoundly affect the shape of geological structures and planetary topography. Accurate descriptions of brittle deformation processes are therefore key to answer fundamental questions such as: How and when does a new fault break? How long can it stay active and under what conditions can tectonic stresses reactivate previously active faults? Which mechanisms promote brittle strain localization and modulate off-fault deformation?

Laboratory experiments have long been used to learn about rock deformation mechanisms in the brittle regime (Paterson & Wong, 2005). The brittle behavior of low-porosity crustal rocks (Figure 1) has some defining characteristics. First and foremost, the differential stress that must be applied to break a rock (the rock's strength) increases with pressure (Byerlee, 1967) (Figure 1a, squares and circles). The stress required to slip on a pre-existing discontinuity is also pressure dependent, and both stresses weakly depend on lithology (Byerlee, 1978). The contrast between these two stresses (intact vs. pre-cut) is typically on the order of hundreds of MPas (Figure 1a). Experiments further reveal a number of phenomena that precede macroscopic failure of a rock sample, such as: a reduction of effective elastic moduli, volume expansion, and acoustic emissions (Figures 1b–1d). Failure is a catastrophic phenomenon that occurs when stresses reach a peak strength which is greater when the imposed strain rate is faster (e.g., Lockner, 1998; Paterson & Wong, 2005). Failure manifests as a transition from distributed to localized strain along macroscopic fractures oriented in a systematic manner with respect to the stress field. It is also well documented that rocks can creep when subjected to a constant stress below their peak strength (e.g., Baud & Meredith, 1997; Brantut et al., 2013; Carter et al., 1981; Heap et al., 2009; Kranz, 1979; Figure 1d). Such brittle creep is typically described as involving three phases: A first phase (primary creep) where strain rate decelerates, a prolonged second phase (secondary creep), during which creep rate remains
# Geochemistry, Geophysics, Geosystems®

## RESEARCH ARTICLE

10.1029/2023GC011229

**Special Section:**

Frontiers in lithospheric dynamics: bridging scales through observations, experiments, and computations

**Key Points:**

* New brittle constitutive law describes the onset of faulting in tectonic simulations

* Model is based on sub-critical growth and interaction of micro-cracks

* Laboratory-derived model parameters can be used to model crustal-scale faulting

**Correspondence to:**

J.-A. Olive,
olive@geologie.ens.fr

**Citation:**

Petit, L., Olive, J.-A., Schubnel, A., Le Pourhiet, L., & Bhat, H. S. (2024). A brittle constitutive law for long-term tectonic modeling based on sub-critical crack growth. *Geochemistry, Geophysics, Geosystems*, 25, e2023GC011229. https://doi.org/10.1029/2023GC011229

Received 7 SEP 2023
Accepted 9 FEB 2024

**Author Contributions:**

**Conceptualization:** Léo Petit, Jean-Arthur Olive, Alexandre Schubnel, Laetitia Le Pourhiet

**Data curation:** Léo Petit

**Formal analysis:** Léo Petit, Jean-Arthur Olive, Alexandre Schubnel, Harsha S. Bhat

**Funding acquisition:** Jean-Arthur Olive, Harsha S. Bhat

**Investigation:** Léo Petit, Jean-Arthur Olive, Alexandre Schubnel, Laetitia Le Pourhiet, Harsha S. Bhat

© 2024 The Authors. Geochemistry, Geophysics, Geosystems published by Wiley Periodicals LLC on behalf of American Geophysical Union. This is an open access article under the terms of the Creative Commons Attribution-NonCommercial-NoDerivs License, which permits use and distribution in any medium, provided the original work is properly cited, the use is non-commercial and no modifications or adaptations are made.

## A Brittle Constitutive Law for Long-Term Tectonic Modeling Based on Sub-Critical Crack Growth

<figure><img src="image_2.png" /></figure>

Léo Petit¹ ID, Jean-Arthur Olive¹ ID, Alexandre Schubnel¹ ID, Laetitia Le Pourhiet² ID, and Harsha S. Bhat¹ ID

¹Laboratoire de Géologie, CNRS - École normale supérieure - PSL University, Paris, France, ²Sorbonne Université, ISTEP, Paris, France

**Abstract** Adequate representations of brittle deformation (fracturing and faulting) are essential ingredients of long-term tectonic simulations. Such models commonly rely on Mohr-Coulomb plasticity coupled with prescribed softening of cohesion and/or friction with accumulated plastic strain. This approach captures fundamental properties of brittle failure, but is overly sensitive to empirical softening parameters that cannot be determined experimentally. Here we design a brittle constitutive law that captures key processes of brittle deformation, and can be straightforwardly implemented in standard geodynamic models. In our Sub-Critically-Altered Maxwell (SCAM) flow law, brittle failure begins with the accumulation of distributed brittle damage, which represents the sub-critical lengthening of tensile micro-cracks prompted by slip on pre-existing shear defects. Damage progressively and permanently weakens the rock's elastic moduli, until cracks catastrophically interact and coalesce up to macroscopic failure. The model's micromechanical parameters can be fully calibrated against rock deformation experiments, alleviating the need for ad-hoc softening parameters. Upon implementing the SCAM flow law in 2-D plane strain simulations of rock deformation experiments, we find that it can produce Coulomb-oriented shear bands which originate as damage bands. SCAM models can also be used to extrapolate rock strength from laboratory to tectonic strain rates, and nuance the use of Byerlee's law as an upper bound on lithosphere stresses. We further show that SCAM models can be upscaled to simulate tectonic deformation of a 10-km thick brittle plate over millions of years. These features make the SCAM rheology a promising tool to further investigate the complexity of brittle behavior across scales.

## 1. Introduction

Tectonic plates tend to be almost rigid and primarily deform within narrow boundary zones. In the upper crust (above ~15 km depth), deformation occurs in the brittle regime through nucleation and growth of fractures and faults, which profoundly affect the shape of geological structures and planetary topography. Accurate descriptions of brittle deformation processes are therefore key to answer fundamental questions such as: How and when does a new fault break? How long can it stay active and under what conditions can tectonic stresses reactivate previously active faults? Which mechanisms promote brittle strain localization and modulate off-fault deformation?

Laboratory experiments have long been used to learn about rock deformation mechanisms in the brittle regime (Paterson & Wong, 2005). The brittle behavior of low-porosity crustal rocks (Figure 1) has some defining characteristics. First and foremost, the differential stress that must be applied to break a rock (the rock's strength) increases with pressure (Byerlee, 1967) (Figure 1a, squares and circles). The stress required to slip on a pre-existing discontinuity is also pressure dependent, and both stresses weakly depend on lithology (Byerlee, 1978). The contrast between these two stresses (intact vs. pre-cut) is typically on the order of hundreds of MPas (Figure 1a). Experiments further reveal a number of phenomena that precede macroscopic failure of a rock sample, such as: a reduction of effective elastic moduli, volume expansion, and acoustic emissions (Figures 1b-1d). Failure is a catastrophic phenomenon that occurs when stresses reach a peak strength which is greater when the imposed strain rate is faster (e.g., Lockner, 1998; Paterson & Wong, 2005). Failure manifests as a transition from distributed to localized strain along macroscopic fractures oriented in a systematic manner with respect to the stress field. It is also well documented that rocks can creep when subjected to a constant stress below their peak strength (e.g., Baud & Meredith, 1997; Brantut et al., 2013; Carter et al., 1981; Heap et al., 2009; Kranz, 1979; Figure 1d). Such brittle creep is typically described as involving three phases: A first phase (primary creep) where strain rate decelerates, a prolonged second phase (secondary creep), during which creep rate remains
<figure><img src="image_3.png" /><figcaption>**Figure 1.** (a) Experimental constraints on brittle rock strength as a function of increasing pressure for intact (squares and circles) and pre-cut (crosses) samples. (b–d) Schematic illustration of axial stress, volumetric strain and acoustic energy versus axial strain, in triaxial experiments performed at a constant strain rate. Point C' marks the onset of dilatancy. Sample cross sections showing the spatial pattern of acoustic emissions are reproduced from Lockner et al. (1991). (e) Typical pattern of axial strain versus time in a brittle creep experiment in which the differential stress $\Delta\sigma$ is imposed and maintained. Colors correspond to values of $\Delta\sigma$ close to (red) or far from (yellow) the rock's peak strength. Here compressive axial strains, strain rates and stresses are plotted as positive numbers for clarity.</figcaption></figure>

nearly steady, and a final stage (tertiary creep) when deformation accelerates until macroscopic failure (Figure 1e, note that the log time representation does not convey the long duration of the secondary phase).

This seemingly complex phenomenology is reasonably well understood as the macroscopic manifestation of the growth and interaction of microcracks that nucleated on pre-existing defects (Tapponnier & Brace, 1976). Crack growth first occurs in a distributed fashion across the sample (Figure 1b). Macroscopic failure then results from the sudden coalescence of interacting microcracks (Figure 1b), whose growth is enabled by differential stress (e.g., Lockner, 1998; Lockner et al., 1991; McBeck et al., 2019). Sample dilatancy points to the tensile nature (mode-I) of some of these cracks, which are susceptible to radiate acoustic energy as they grow (Figure 1d). The time and strain rate dependence of these phenomena further suggests that the speed of crack propagation in the
<figure><img src="image_3.png" /><figcaption>**Figure 1.** (a) Experimental constraints on brittle rock strength as a function of increasing pressure for intact (squares and circles) and pre-cut (crosses) samples. (b–d) Schematic illustration of axial stress, volumetric strain and acoustic energy versus axial strain, in triaxial experiments performed at a constant strain rate. Point C' marks the onset of dilatancy. Sample cross sections showing the spatial pattern of acoustic emissions are reproduced from Lockner et al. (1991). (e) Typical pattern of axial strain versus time in a brittle creep experiment in which the differential stress $\Delta\sigma$ is imposed and maintained. Colors correspond to values of $\Delta\sigma$ close to (red) or far from (yellow) the rock's peak strength. Here compressive axial strains, strain rates and stresses are plotted as positive numbers for clarity.</figcaption></figure>

nearly steady, and a final stage (tertiary creep) when deformation accelerates until macroscopic failure (Figure 1e, note that the log time representation does not convey the long duration of the secondary phase).

This seemingly complex phenomenology is reasonably well understood as the macroscopic manifestation of the growth and interaction of microcracks that nucleated on pre-existing defects (Tapponnier & Brace, 1976). Crack growth first occurs in a distributed fashion across the sample (Figure 1b). Macroscopic failure then results from the sudden coalescence of interacting microcracks (Figure 1b), whose growth is enabled by differential stress (e.g., Lockner, 1998; Lockner et al., 1991; McBeck et al., 2019). Sample dilatancy points to the tensile nature (mode-I) of some of these cracks, which are susceptible to radiate acoustic energy as they grow (Figure 1d). The time and strain rate dependence of these phenomena further suggests that the speed of crack propagation in the
bulk rock depends on the forces acting at crack tips, which is typical of sub-critical crack growth processes. The main underlying mechanism in the brittle regime is known as stress corrosion (Atkinson, 1984). It refers to reactions occurring between a chemically active fluid and the strained atomic bonds at the tip of microcracks, which induce stress-dependent kinetics of bond breaking (Eppes & Keanini, 2017). While other mechanisms such as pressure-solution (e.g., Gratier et al., 2013), can also contribute to rate-dependent deformation as pressure increases, sub-critical crack growth has been identified as a key contributor to the strain rate (i.e., time-) dependent behavior of brittle rocks in the brittle regime that is particularly well highlighted by brittle creep experiments (Brantut et al., 2013).

Though the phenomenology of brittle failure was well known long before geodynamicists harnessed the power of microprocessors, most tectonic simulations currently rely on a simplified treatment that consists in capping stresses at a rate-independent Mohr-Coulomb yield envelope (e.g., Gerya, 2010; Poliakov & Buck, 1998). This has the advantage of being numerically efficient, adequately capturing the pressure-dependent frictional strength of pre-cut rocks, and spontaneously localizing plastic strain through the bifurcative properties of the Mohr-Coulomb plastic flow rule (e.g., Kaus, 2010; Lemiale et al., 2008; Rudnicki & Rice, 1975; Vermeer & De Borst, 1984). In this framework, strain localization is typically accompanied by a rotation of the principal stresses inside the incipient shear band, which leads to a reduction of the remote stresses (Le Pourhiet, 2013). By itself, this rotation-induced “structural” softening does not account for the 100 s of MPas that separate the strength of intact rocks from their residual strength once faulted (Figure 1a). An approach commonly used to promote sustained strain localization in tectonic simulations (Figure 2) is to weaken the material friction μ and cohesion C, from {μmax, Cmax} to {μmin, Cmin} over a certain amount of non-recoverable (plastic) strain Δe<sub>II</sub><sup>p</sup> (e.g., Lavier et al., 2000; Poliakov & Buck, 1998; Figure 2a). This amounts to enforcing a contrast between intact and broken rocks reminiscent of the strength contrast observed experimentally.

Strain-weakened Mohr-Coulomb plasticity however presents several drawbacks. This parameterization typically ignores the strain rate dependence of rocks' intact strength, and relies on a single value of intact friction and cohesion to determine the intact yield strength. Further, the critical plastic strain Δe<sub>II</sub><sup>p</sup> is meant to represent a wide range of possible weakening mechanisms, and is therefore not easily quantified through laboratory experiments. These limitations can be problematic since the choice of weakening parameters can have major consequences on the outcome of a tectonic simulation. Lavier et al. (2000) for example, pointed out the spectacular effect of Δe<sub>II</sub><sup>p</sup> on tectonic styles produced in a rifting simulation (Figures 2b vs. 2c). While some recent studies have investigated the effects of various weakening parameterizations (e.g., Duretz et al. (2021); Meyer et al. (2017); Naliboff

<figure><img src="image_4.png" /><figcaption>**Figure 2.** (a) Schematic representation of an elastic-plastic rheology with strain weakening, under constant applied strain rate. The difference between initial ($\sigma_y^{(MAX)}$) and final ($\sigma_y^{(MIN)}$) yield stresses are caused by a prescribed decrease of frictional properties $\mu$ and $C$ over a specified amount of plastic strain $\Delta e_{II}^p$. (b–c) Example simulation of extension in a 10 km-thick elastic-plastic upper crust overlying an inviscid medium (Lavier et al., 2000). The only difference between the two panels is the choice of $\Delta e_{II}^p$, which is small in b, producing a large-offset normal fault (detachment), and large in c, producing two conjugate faults that outline a graben structure. In this example the difference between $\sigma_y^{(MAX)}$ and $\sigma_y^{(MIN)}$ is caused by a drop in material cohesion while friction is kept constant.</figcaption></figure>
bulk rock depends on the forces acting at crack tips, which is typical of sub-critical crack growth processes. The main underlying mechanism in the brittle regime is known as stress corrosion (Atkinson, 1984). It refers to reactions occurring between a chemically active fluid and the strained atomic bonds at the tip of microcracks, which induce stress-dependent kinetics of bond breaking (Eppes & Keanini, 2017). While other mechanisms such as pressure-solution (e.g., Gratier et al., 2013), can also contribute to rate-dependent deformation as pressure increases, sub-critical crack growth has been identified as a key contributor to the strain rate (i.e., time-) dependent behavior of brittle rocks in the brittle regime that is particularly well highlighted by brittle creep experiments (Brantut et al., 2013).

Though the phenomenology of brittle failure was well known long before geodynamicists harnessed the power of microprocessors, most tectonic simulations currently rely on a simplified treatment that consists in capping stresses at a rate-independent Mohr-Coulomb yield envelope (e.g., Gerya, 2010; Poliakov & Buck, 1998). This has the advantage of being numerically efficient, adequately capturing the pressure-dependent frictional strength of pre-cut rocks, and spontaneously localizing plastic strain through the bifurcative properties of the Mohr-Coulomb plastic flow rule (e.g., Kaus, 2010; Lemiale et al., 2008; Rudnicki & Rice, 1975; Vermeer & De Borst, 1984). In this framework, strain localization is typically accompanied by a rotation of the principal stresses inside the incipient shear band, which leads to a reduction of the remote stresses (Le Pourhiet, 2013). By itself, this rotation-induced “structural” softening does not account for the 100 s of MPas that separate the strength of intact rocks from their residual strength once faulted (Figure 1a). An approach commonly used to promote sustained strain localization in tectonic simulations (Figure 2) is to weaken the material friction μ and cohesion C, from {μmax, Cmax} to {μmin, Cmin} over a certain amount of non-recoverable (plastic) strain Δe<sub>II</sub><sup>p</sup> (e.g., Lavier et al., 2000; Poliakov & Buck, 1998; Figure 2a). This amounts to enforcing a contrast between intact and broken rocks reminiscent of the strength contrast observed experimentally.

Strain-weakened Mohr-Coulomb plasticity however presents several drawbacks. This parameterization typically ignores the strain rate dependence of rocks' intact strength, and relies on a single value of intact friction and cohesion to determine the intact yield strength. Further, the critical plastic strain Δe<sub>II</sub><sup>p</sup> is meant to represent a wide range of possible weakening mechanisms, and is therefore not easily quantified through laboratory experiments. These limitations can be problematic since the choice of weakening parameters can have major consequences on the outcome of a tectonic simulation. Lavier et al. (2000) for example, pointed out the spectacular effect of Δe<sub>II</sub><sup>p</sup> on tectonic styles produced in a rifting simulation (Figures 2b vs. 2c). While some recent studies have investigated the effects of various weakening parameterizations (e.g., Duretz et al. (2021); Meyer et al. (2017); Naliboff

<figure><img src="image_4.png" /><figcaption>**Figure 2.** (a) Schematic representation of an elastic-plastic rheology with strain weakening, under constant applied strain rate. The difference between initial ($\sigma_y^{(MAX)}$) and final ($\sigma_y^{(MIN)}$) yield stresses are caused by a prescribed decrease of frictional properties $\mu$ and $C$ over a specified amount of plastic strain $\Delta e_{II}^p$. (b–c) Example simulation of extension in a 10 km-thick elastic-plastic upper crust overlying an inviscid medium (Lavier et al., 2000). The only difference between the two panels is the choice of $\Delta e_{II}^p$, which is small in b, producing a large-offset normal fault (detachment), and large in c, producing two conjugate faults that outline a graben structure. In this example the difference between $\sigma_y^{(MAX)}$ and $\sigma_y^{(MIN)}$ is caused by a drop in material cohesion while friction is kept constant.</figcaption></figure>
et al. (2020); Pan et al. (2023)), it remains common practice to rely on ad-hoc softening rules in geodynamic simulations without assessing their impact on model behavior.

One path toward remedying this issue is to improve the way geodynamic simulations parameterize the transition from intact to broken rock, in a manner that allows more direct comparison with experimental data and can be interpreted in terms of underlying deformation mechanisms. An adequate parameterization of progressive brittle failure should indeed account for standard observations such as the pre-peak reduction in elastic moduli, the evolving spatial pattern of acoustic emissions, or sample dilatancy which ceases upon failure (Figures 1b–1d). It should also account for the strain rate dependence of brittle yielding and the occurrence of brittle creep. Finally, it should include a representation of the ever-evolving internal state of the rock to include a memory of past deformation events. A promising alternative is to turn to models that describe brittle yielding as the accumulation of damage which ultimately leads to macroscopic failure.

A first family of such models are Continuum Damage Mechanics models. They treat failure as a progressive phenomenon indexed on the alteration of a rock's internal state (damage), and can produce strain rate-dependent brittle strengths, as well as pre-peak softening. Some are built on thermodynamic descriptions of energy dissipation during inelastic deformation (e.g., Hamiel et al., 2004; Karrech et al., 2011; Lyakhovsky et al., 1997), others simply index damage growth on excess stresses above a yield stress, and strain (e.g., Manaker et al., 2006). They do not assume a specific microstructure, which makes them flexible but also not directly interpretable in terms of deformation processes.

In that regard, micromechanics-based models have been particularly successful at capturing the broad range of behaviors associated with brittle deformation (Paterson & Wong, 2005). In this family of damage models, assumptions about the distribution and geometry of pre-existing defects in the material allow the analytical determination of stress concentrations around them, using linear elastic fracture mechanics. Motion along defects cause the stress intensity factors (i.e., a measure of the stress state at the edge of discontinuities) at their tips to increase up to the fracture toughness of the rock, allowing tensile crack propagation. Drivers of such stress heterogeneities can be planar flaws such as grain boundaries, pre-existing microcracks (e.g., Ashby & Sammis, 1990; Kachanov, 1982a, 1982c; Nemat-Nasser & Horii, 1982), pores (Sammis & Ashby, 1986), moduli contrasts across grains in contact (Dey & Wang, 1981), or can even be left undetermined (e.g., Costin, 1985). Tensile cracks, in turn, alter the effective elastic properties of the rock as they lengthen, in an anisotropic fashion (Budiansky & O'connell, 1976; Deshpande & Evans, 2008; Kachanov, 1993; Walsh, 1965a, 1965b). This framework has been used to model high strain rate deformation (e.g., during seismic rupture, Bhat et al., 2012; Thomas et al., 2017) assuming critical fracture propagation, as well as slow deformation assuming sub-critical crack growth (Kachanov, 1982b). The latter class of models has also been used to describe brittle creep, assuming pre-existing planar defects (Brantut et al., 2012), successfully accounting for the multi-phased dynamics of brittle creep (Figure 1e).

One drawback of this approach is its computational cost, because it requires to accurately resolve the kinetics of fracture lengthening, which crack interactions ultimately render unstable close to macroscopic failure. This may explain why it has not yet been implemented in long-term, large scale tectonic simulations, even though the processes it describes are clearly central to the initiation and evolution of crustal faults. By representing specific deformation mechanisms that can be studied in the laboratory, these models can indeed be calibrated against experiments and need not resort to ad-hoc macroscopic parameters (e.g., Bhat et al., 2011; Brantut et al., 2012; Costin, 1983, 1985).

As a first step in this direction, this study aims at constructing a constitutive brittle rheology rooted in the subcritical growth of microcracks from pre-existing rock defects. We seek a formulation that (a) captures the essence of brittle rock behavior at the expense of a few simplifications, (b) has a straightforward micromechanical interpretation, (c) can be calibrated against experimental data, and (d) is usable in standard 2-D plane strain numerical geodynamic models. We propose such a constitutive law in Section 2, and describe its fundamental behavior in terms of stress-strain curves in Section 3. This allows us to calibrate its parameters using experimental data from both constant strain rate and brittle creep tests. We then implement our constitutive law in 2-D plane strain numerical simulations that reproduce experimental conditions (Section 4), and discuss the model's key features in Section 5. Finally, we implement our constitutive law in a crustal-scale tectonic simulator and compare it to the standard elasto-plastic approach (Section 6).
et al. (2020); Pan et al. (2023)), it remains common practice to rely on ad-hoc softening rules in geodynamic simulations without assessing their impact on model behavior.

One path toward remedying this issue is to improve the way geodynamic simulations parameterize the transition from intact to broken rock, in a manner that allows more direct comparison with experimental data and can be interpreted in terms of underlying deformation mechanisms. An adequate parameterization of progressive brittle failure should indeed account for standard observations such as the pre-peak reduction in elastic moduli, the evolving spatial pattern of acoustic emissions, or sample dilatancy which ceases upon failure (Figures 1b–1d). It should also account for the strain rate dependence of brittle yielding and the occurrence of brittle creep. Finally, it should include a representation of the ever-evolving internal state of the rock to include a memory of past deformation events. A promising alternative is to turn to models that describe brittle yielding as the accumulation of damage which ultimately leads to macroscopic failure.

A first family of such models are Continuum Damage Mechanics models. They treat failure as a progressive phenomenon indexed on the alteration of a rock's internal state (damage), and can produce strain rate-dependent brittle strengths, as well as pre-peak softening. Some are built on thermodynamic descriptions of energy dissipation during inelastic deformation (e.g., Hamiel et al., 2004; Karrech et al., 2011; Lyakhovsky et al., 1997), others simply index damage growth on excess stresses above a yield stress, and strain (e.g., Manaker et al., 2006). They do not assume a specific microstructure, which makes them flexible but also not directly interpretable in terms of deformation processes.

In that regard, micromechanics-based models have been particularly successful at capturing the broad range of behaviors associated with brittle deformation (Paterson & Wong, 2005). In this family of damage models, assumptions about the distribution and geometry of pre-existing defects in the material allow the analytical determination of stress concentrations around them, using linear elastic fracture mechanics. Motion along defects cause the stress intensity factors (i.e., a measure of the stress state at the edge of discontinuities) at their tips to increase up to the fracture toughness of the rock, allowing tensile crack propagation. Drivers of such stress heterogeneities can be planar flaws such as grain boundaries, pre-existing microcracks (e.g., Ashby & Sammis, 1990; Kachanov, 1982a, 1982c; Nemat-Nasser & Horii, 1982), pores (Sammis & Ashby, 1986), moduli contrasts across grains in contact (Dey & Wang, 1981), or can even be left undetermined (e.g., Costin, 1985). Tensile cracks, in turn, alter the effective elastic properties of the rock as they lengthen, in an anisotropic fashion (Budiansky & O'connell, 1976; Deshpande & Evans, 2008; Kachanov, 1993; Walsh, 1965a, 1965b). This framework has been used to model high strain rate deformation (e.g., during seismic rupture, Bhat et al., 2012; Thomas et al., 2017) assuming critical fracture propagation, as well as slow deformation assuming sub-critical crack growth (Kachanov, 1982b). The latter class of models has also been used to describe brittle creep, assuming pre-existing planar defects (Brantut et al., 2012), successfully accounting for the multi-phased dynamics of brittle creep (Figure 1e).

One drawback of this approach is its computational cost, because it requires to accurately resolve the kinetics of fracture lengthening, which crack interactions ultimately render unstable close to macroscopic failure. This may explain why it has not yet been implemented in long-term, large scale tectonic simulations, even though the processes it describes are clearly central to the initiation and evolution of crustal faults. By representing specific deformation mechanisms that can be studied in the laboratory, these models can indeed be calibrated against experiments and need not resort to ad-hoc macroscopic parameters (e.g., Bhat et al., 2011; Brantut et al., 2012; Costin, 1983, 1985).

As a first step in this direction, this study aims at constructing a constitutive brittle rheology rooted in the subcritical growth of microcracks from pre-existing rock defects. We seek a formulation that (a) captures the essence of brittle rock behavior at the expense of a few simplifications, (b) has a straightforward micromechanical interpretation, (c) can be calibrated against experimental data, and (d) is usable in standard 2-D plane strain numerical geodynamic models. We propose such a constitutive law in Section 2, and describe its fundamental behavior in terms of stress-strain curves in Section 3. This allows us to calibrate its parameters using experimental data from both constant strain rate and brittle creep tests. We then implement our constitutive law in 2-D plane strain numerical simulations that reproduce experimental conditions (Section 4), and discuss the model's key features in Section 5. Finally, we implement our constitutive law in a crustal-scale tectonic simulator and compare it to the standard elasto-plastic approach (Section 6).
## 2. A Sub-Critically Altered Maxwell (SCAM) Constitutive Law for Brittle Deformation

### 2.1. Generic Stress-Strain Relation

Our constitutive model builds upon an isotropic, incompressible elastic stress-strain relationship:

$$e_{ij} = \frac{s_{ij}}{2G}, \qquad (1)$$

linking the deviatoric strain and stress tensors, $e_{ij}$ and $s_{ij}$, through shear modulus $G$. In the following, we adopt the convention of summed repeated indices. Our fundamental assumption is that the shear modulus is altered as a function of the internal state of the material, which leads to path-dependent behavior. Specifically, we assume that $G$ decreases as a function of a scalar state variable $D$, a measure of rock damage, to be defined in Section 2.2:

$$G = G_0 f(D). \qquad (2)$$

In Equation 2, $G_0$ denotes the shear modulus of the material in its least damaged state, and $f(D)$ a decreasing scalar function of $D$, hereafter referred to as “weakening function,” satisfying $f(D) \in [0,1]$ and $f'(D) < 0$. The incompressible elastic relationship Equation 1 can be recast as a damaged-elastic constitutive law

$$e_{ij} = \frac{s_{ij}}{2G_0 f(D)}, \qquad (3)$$

which takes the form of a Maxwell visco-elastic constitutive law upon time differentiation:

$$\dot{e}_{ij} = \frac{\dot{s}_{ij}}{2G_0 f(D)} + \frac{s_{ij}}{2\eta_D}. \qquad (4)$$

In Equation 4 $\eta_D$ is a viscosity associated with damage growth:

$$\eta_D = \frac{f^2(D) G_0}{|f'(D)| \dot{D}}. \qquad (5)$$

The damage state variable is related to the lengthening of mode-I microfractures, an intrinsically dilatant process. Throughout this study, inelastic dilatancy is neglected in favor of a purely deviatoric description of the damaged rheology, focusing on the role of microcracking on shear modulus alteration, and on fault nucleation. A strategy to account for damage-induced dilatancy within the SCAM framework will nonetheless be outlined in Section 7. In the following, we detail the micromechanical interpretation of the damage variable, the model governing its growth rate, as well as the weakening function.

### 2.2. Micromechanical Representation of Rock Damage

Our goal is to model the accumulation of damage in the upper crust, which is primarily composed of low-porosity (<1%) magmatic and metamorphic silicate rocks. These units lie in an overall compressive stress state, with pressures up to hundreds of MPas. Yet, distributed brittle deformation typically involves the opening of mode-I microcracks (Figure 1), which is made possible by stress concentrations around defects or grain boundaries. To describe these processes, we adopt the damage framework developed by Ashby and Sammis (1990) which has been used successfully to predict the brittle strength of several rocks (e.g., Bhat et al., 2011; Baud et al., 2000; X. Wu et al., 2000) at low confining pressure, and the dynamics of fracturing during seismic ruptures (Bhat et al., 2012; Thomas et al., 2017). This model considers the growth of tensile “wing”-cracks from the tips of penny-shaped shear defects distributed within the rock (Figures 3a–3d).

The damage variable represents the relative volume occupied by cracks as the wings lengthen in the direction of the most compressive stress (Figure 3d). It is defined as
## 2. A Sub-Critically Altered Maxwell (SCAM) Constitutive Law for Brittle Deformation

### 2.1. Generic Stress-Strain Relation

Our constitutive model builds upon an isotropic, incompressible elastic stress-strain relationship:

$$e_{ij} = \frac{s_{ij}}{2G}, \qquad (1)$$

linking the deviatoric strain and stress tensors, $e_{ij}$ and $s_{ij}$, through shear modulus $G$. In the following, we adopt the convention of summed repeated indices. Our fundamental assumption is that the shear modulus is altered as a function of the internal state of the material, which leads to path-dependent behavior. Specifically, we assume that $G$ decreases as a function of a scalar state variable $D$, a measure of rock damage, to be defined in Section 2.2:

$$G = G_0 f(D). \qquad (2)$$

In Equation 2, $G_0$ denotes the shear modulus of the material in its least damaged state, and $f(D)$ a decreasing scalar function of $D$, hereafter referred to as “weakening function,” satisfying $f(D) \in [0,1]$ and $f'(D) < 0$. The incompressible elastic relationship Equation 1 can be recast as a damaged-elastic constitutive law

$$e_{ij} = \frac{s_{ij}}{2G_0 f(D)}, \qquad (3)$$

which takes the form of a Maxwell visco-elastic constitutive law upon time differentiation:

$$\dot{e}_{ij} = \frac{\dot{s}_{ij}}{2G_0 f(D)} + \frac{s_{ij}}{2\eta_D}. \qquad (4)$$

In Equation 4 $\eta_D$ is a viscosity associated with damage growth:

$$\eta_D = \frac{f^2(D) G_0}{|f'(D)| \dot{D}}. \qquad (5)$$

The damage state variable is related to the lengthening of mode-I microfractures, an intrinsically dilatant process. Throughout this study, inelastic dilatancy is neglected in favor of a purely deviatoric description of the damaged rheology, focusing on the role of microcracking on shear modulus alteration, and on fault nucleation. A strategy to account for damage-induced dilatancy within the SCAM framework will nonetheless be outlined in Section 7. In the following, we detail the micromechanical interpretation of the damage variable, the model governing its growth rate, as well as the weakening function.

### 2.2. Micromechanical Representation of Rock Damage

Our goal is to model the accumulation of damage in the upper crust, which is primarily composed of low-porosity (<1%) magmatic and metamorphic silicate rocks. These units lie in an overall compressive stress state, with pressures up to hundreds of MPas. Yet, distributed brittle deformation typically involves the opening of mode-I microcracks (Figure 1), which is made possible by stress concentrations around defects or grain boundaries. To describe these processes, we adopt the damage framework developed by Ashby and Sammis (1990) which has been used successfully to predict the brittle strength of several rocks (e.g., Bhat et al., 2011; Baud et al., 2000; X. Wu et al., 2000) at low confining pressure, and the dynamics of fracturing during seismic ruptures (Bhat et al., 2012; Thomas et al., 2017). This model considers the growth of tensile “wing”-cracks from the tips of penny-shaped shear defects distributed within the rock (Figures 3a–3d).

The damage variable represents the relative volume occupied by cracks as the wings lengthen in the direction of the most compressive stress (Figure 3d). It is defined as
<figure><img src="image_5.png" /><figcaption>**Figure 3.** (a–c) Cross-sections in the ($\sigma_1$, $\sigma_3$) plane showing the array of cracks growing in a rock, in the model of Ashby and Sammis (1990). These cracks initially grow from a shear defect of radius $a$, oriented at an angle $\psi$ to $\sigma_1$. At its tips, “wing” tensile cracks of length $l$ (d) may preexist. Under differential stresses unable to overcome the frictional resistance along shear defects (a), the stress concentration ($K_I$) at their edges is negative (e) and the material deforms elastically. Once stresses overcome the frictional resistance, the sliding shear defects exert a wedging force $F_w$ increasing $K_I$ (f). Positive $K_I$ promotes stable tensile cracks growth (b). As cracks lengthen, they begin to interact (c). At this point damage growth enters an unstable, accelerating regime (g). (h) Alteration of the shear modulus as a function of damage.</figcaption></figure>

$$D = \frac{4}{3}\pi N_v(\alpha a + l)^3, \qquad (6)$$

where $N_v$ is the number of shear defects per unit volume, $a$ and $l$ are the radius of the shear defects and the length of the wing cracks, respectively. $\alpha = \cos \psi$ is the cosine of the angle $\psi$ between the shear defects and the most compressive stress (Figure 3d). The least damaged state ($D = D_0$) corresponds to the state of a rock containing only shear defects, where no wing crack has nucleated (i.e., $l = 0$). The most damaged state occurs at $D = 1$ when the volume of the spheres enclosing each wing crack has grown to match the characteristic volume defined by the spacing of defects ($V_c = 1/N_v$). This upper bound is the result of the formulation of the interaction between cracks, detailed in Section 2.3. It represents a stage at which coalescence of cracks becomes unavoidable. For simplicity, we only consider shear defects with normal vectors lying in the plane of the two extreme principal stresses $\sigma_1$ and $\sigma_3$. This allows us to index their activity on a 2-D Mohr-Coulomb yield criterion.

## 2.3. Damage Growth

We assume that under low strain rates and on long time scales, wing cracks lengthen in a sub-critical manner, that is, with stress intensity factors ($K_I$) lower than the fracture toughness ($K_{IC}$) of the material (Atkinson, 1984). To capture this process in our constitutive law, we adopt the stress corrosion law introduced by Charles (1958), which has proven successful at explaining experimental data (e.g., Atkinson, 1984; Brantut et al., 2012; Deshpande & Evans, 2008; Kachanov, 1982c). Specifically, the crack growth rate writes

$$\dot{i} = \dot{i}_0 \left( \frac{K_I}{K_{IC}} \right)^n, \qquad (7)$$
<figure><img src="image_5.png" /><figcaption>**Figure 3.** (a–c) Cross-sections in the ($\sigma_1$, $\sigma_3$) plane showing the array of cracks growing in a rock, in the model of Ashby and Sammis (1990). These cracks initially grow from a shear defect of radius $a$, oriented at an angle $\psi$ to $\sigma_1$. At its tips, “wing” tensile cracks of length $l$ (d) may preexist. Under differential stresses unable to overcome the frictional resistance along shear defects (a), the stress concentration ($K_I$) at their edges is negative (e) and the material deforms elastically. Once stresses overcome the frictional resistance, the sliding shear defects exert a wedging force $F_w$ increasing $K_I$ (f). Positive $K_I$ promotes stable tensile cracks growth (b). As cracks lengthen, they begin to interact (c). At this point damage growth enters an unstable, accelerating regime (g). (h) Alteration of the shear modulus as a function of damage.</figcaption></figure>

$$D = \frac{4}{3}\pi N_v(\alpha a + l)^3, \qquad (6)$$

where $N_v$ is the number of shear defects per unit volume, $a$ and $l$ are the radius of the shear defects and the length of the wing cracks, respectively. $\alpha = \cos \psi$ is the cosine of the angle $\psi$ between the shear defects and the most compressive stress (Figure 3d). The least damaged state ($D = D_0$) corresponds to the state of a rock containing only shear defects, where no wing crack has nucleated (i.e., $l = 0$). The most damaged state occurs at $D = 1$ when the volume of the spheres enclosing each wing crack has grown to match the characteristic volume defined by the spacing of defects ($V_c = 1/N_v$). This upper bound is the result of the formulation of the interaction between cracks, detailed in Section 2.3. It represents a stage at which coalescence of cracks becomes unavoidable. For simplicity, we only consider shear defects with normal vectors lying in the plane of the two extreme principal stresses $\sigma_1$ and $\sigma_3$. This allows us to index their activity on a 2-D Mohr-Coulomb yield criterion.

## 2.3. Damage Growth

We assume that under low strain rates and on long time scales, wing cracks lengthen in a sub-critical manner, that is, with stress intensity factors ($K_I$) lower than the fracture toughness ($K_{IC}$) of the material (Atkinson, 1984). To capture this process in our constitutive law, we adopt the stress corrosion law introduced by Charles (1958), which has proven successful at explaining experimental data (e.g., Atkinson, 1984; Brantut et al., 2012; Deshpande & Evans, 2008; Kachanov, 1982c). Specifically, the crack growth rate writes

$$\dot{i} = \dot{i}_0 \left( \frac{K_I}{K_{IC}} \right)^n, \qquad (7)$$
where $\dot{l}_0$ is a characteristic crack growth rate and n the Charles law exponent. The damage growth rate is then retrieved from the wing-crack tip speed

$$ \dot{D} = \frac{\partial D}{\partial l} l = \frac{3D_0^2 D_0^{\frac{1}{2}} \dot{l}_0}{\alpha a (K_I K_{IC})^n} \quad (8) $$

This equation applies only when $K_I > 0$, otherwise $\dot{D} = 0$. Using Equation 7 requires an expression for $K_I$, the stress intensity factor at the tip of the wing cracks. Following Ashby and Sammis (1990), we assemble $K_I$ as the sum of three terms:

$$ K_I = K_I^{(w)} + K_I^{(\sigma_3)} + K_I^{(i)} \quad (9) $$

The first term ($K_I^{(w)}$) represents the stress concentration due to frictional slip on the shear defects wedging open the wing cracks. Following Tada et al. (1973), it can be expressed as the action of a tensile wedging force $F_w$ at the center of an equivalent penny-shaped crack. The radius of this circular crack is that of the sphere enclosing one entire wing crack (shear defect + tensile wings, Figure 3d). However, instead of writing it $l + \alpha a$, as in Equation 6, we write it $l + \beta a$, where $\beta$ is a regularization factor. This approach was adopted by Ashby and Sammis (1990) to ensure that in the absence of wing cracks ($l = 0, D = D_0$), $K_I$ matches the stress intensity factor at the tip of shear defects as derived by Ashby and Hallam (1986). This yields $\beta = 1/\pi$ (Bhat et al., 2011) and the following expression for $K_I^{(w)}$:

$$ K_I^{(w)} = \frac{F_w}{[\pi(l + \beta a)]^{3/2}} \quad (10) $$

The wedging force relates to the excess shear stress acting on the defects (of area $\pi a^2$) relative to their frictional resistance. Following Ashby and Sammis (1990), we write:

$$ F_w = (\sigma_3 A_3 - \sigma_1 A_1) a^2. \quad (11) $$

$A_1$ and $A_3$ are constants that depend on the friction and orientation of the shear defects. In the following, we assume $\psi = 45^\circ$ as Ashby and Hallam (1986) showed that this orientation maximizes the wedging force over a wide range of wing crack lengths. This yields:

$$ A_1 = \pi \sqrt{\frac{\beta}{3} [\sqrt{1 + \mu^2} - \mu]} \quad (12) $$

$$ A_3 = A_1 \left[ \frac{\sqrt{1 + \mu^2 + \mu}}{\sqrt{1 + \mu^2} - \mu} \right]. \quad (13) $$

Overall, $K_I^{(w)}$ strongly depends on the differential stress $\Delta\sigma = \sigma_3 - \sigma_1$ that develops in the rock, as it allows frictional slip on the defects and wedging of the wings. By contrast, the second term in Equation 9 ($K_I^{(\sigma_3)}$) represents remote wing-normal compression $\sigma_3$ acting to close tensile cracks. Bhat et al. (2011) estimated it based on results from Tada et al. (1973) as:

$$ K_I^{(\sigma_3)} = \frac{2}{\pi} (\sigma_3) \sqrt{\pi l}. \quad (14) $$

Finally, the third term ($K_I^{(i)}$) serves to describe the interaction between cracks as they lengthen, and is a core feature of this micromechanical model. Ashby and Sammis (1990) required that the wedging forces applied to cracks be compensated by an internal stress ($\sigma_i$ in Figure 3) to satisfy mechanical equilibrium. The internal stress is applied on an effective area perpendicular to $\sigma_3$ that separates neighboring cracks ($A_b$). The sum of this area with
where $\dot{l}_0$ is a characteristic crack growth rate and $n$ the Charles law exponent. The damage growth rate is then retrieved from the wing-crack tip speed

$$ \dot{D} = \frac{\partial D}{\partial l} l = \frac{3D_0^2 D_0^{\frac{1}{2}} \dot{l}_0}{\alpha a K_{IC}} \left( \frac{K_I}{K_{IC}} \right)^n \quad (8) $$

This equation applies only when $K_I > 0$, otherwise $\dot{D} = 0$. Using Equation 7 requires an expression for $K_I$, the stress intensity factor at the tip of the wing cracks. Following Ashby and Sammis (1990), we assemble $K_I$ as the sum of three terms:

$$ K_I = K_I^{(w)} + K_I^{(\sigma_3)} + K_I^{(i)} \quad (9) $$

The first term ($K_I^{(w)}$) represents the stress concentration due to frictional slip on the shear defects wedging open the wing cracks. Following Tada et al. (1973), it can be expressed as the action of a tensile wedging force $F_w$ at the center of an equivalent penny-shaped crack. The radius of this circular crack is that of the sphere enclosing one entire wing crack (shear defect + tensile wings, Figure 3d). However, instead of writing it $l + \alpha a$, as in Equation 6, we write it $l + \beta a$, where $\beta$ is a regularization factor. This approach was adopted by Ashby and Sammis (1990) to ensure that in the absence of wing cracks ($l = 0, D = D_0$), $K_I$ matches the stress intensity factor at the tip of shear defects as derived by Ashby and Hallam (1986). This yields $\beta = 1/\pi$ (Bhat et al., 2011) and the following expression for $K_I^{(w)}$:

$$ K_I^{(w)} = \frac{F_w}{[\pi(l + \beta a)]^{3/2}} \quad (10) $$

The wedging force relates to the excess shear stress acting on the defects (of area $\pi a^2$) relative to their frictional resistance. Following Ashby and Sammis (1990), we write:

$$ F_w = (\sigma_3 A_3 - \sigma_1 A_1) a^2. \quad (11) $$

$A_1$ and $A_3$ are constants that depend on the friction and orientation of the shear defects. In the following, we assume $\psi = 45^\circ$ as Ashby and Hallam (1986) showed that this orientation maximizes the wedging force over a wide range of wing crack lengths. This yields:

$$ A_1 = \pi \sqrt{\frac{\beta}{3} [\sqrt{1 + \mu^2} - \mu]} \quad (12) $$

$$ A_3 = A_1 \left[ \frac{\sqrt{1 + \mu^2 + \mu}}{\sqrt{1 + \mu^2} - \mu} \right]. \quad (13) $$

Overall, $K_I^{(w)}$ strongly depends on the differential stress $\Delta\sigma = \sigma_3 - \sigma_1$ that develops in the rock, as it allows frictional slip on the defects and wedging of the wings. By contrast, the second term in Equation 9 ($K_I^{(\sigma_3)}$) represents remote wing-normal compression $\sigma_3$ acting to close tensile cracks. Bhat et al. (2011) estimated it based on results from Tada et al. (1973) as:

$$ K_I^{(\sigma_3)} = \frac{2}{\pi} (\sigma_3) \sqrt{\pi l}. \quad (14) $$

Finally, the third term ($K_I^{(i)}$) serves to describe the interaction between cracks as they lengthen, and is a core feature of this micromechanical model. Ashby and Sammis (1990) required that the wedging forces applied to cracks be compensated by an internal stress ($\sigma_i$ in Figure 3) to satisfy mechanical equilibrium. The internal stress is applied on an effective area perpendicular to $\sigma_3$ that separates neighboring cracks ($A_b$). The sum of this area with
the characteristic area of each wing crack ($\pi(l + \alpha a)^2$) amounts to the area $A_c$ that is obtained by projecting the spherical volume $V_c = 1/N_v$ along $\sigma_3$. Therefore, $A_b = A_c - \pi(l + \alpha a)^2$, with

$$A_c = \pi^{1/3} \left( \frac{3}{4N_v} \right)^{2/3}. \qquad (15)$$

This leads to the following expression for the internal stress acting in the direction of least compression, $\sigma_3^i$:

$$\sigma_3^i = \frac{F_w}{A_b}. \qquad (16)$$

Internal stress $\sigma_3^i$ increases dramatically as wings lengthen (D approaches 1) and the areas between fractures ($A_b$) shrink. This is when crack interactions become dominant. $K_I^{(i)}$ is readily obtained from $\sigma_3^i$ by analogy with Equation 14:

$$K_I^{(i)} = \frac{2}{\pi}(\sigma_3^i)\sqrt{\pi l}. \qquad (17)$$

The full expression of $K_I$ then reads

$$K_I = \frac{F_w}{[\pi(l + \beta a)]^{3/2}} + \frac{2}{\pi}(\sigma_3 + \sigma_3^i)\sqrt{\pi l}. \qquad (18)$$

It can be recast as a function of damage rather than crack length, following Bhat et al. (2011), yielding

$$K_I = \sqrt{\pi a}[(\sigma_3 A_3 - \sigma_1 A_1)(c_1 + c_2) + \sigma_3 c_3], \qquad (19)$$

where $c_1$, $c_2$, and $c_3$ are functions of the damage state that write

$$c_1 = \frac{1}{\pi^2 \alpha^{3/2} [(D/D_0)^{1/3} - 1 + \beta/\alpha]^{3/2}} \qquad (20)$$

$$c_2 = \frac{2}{\pi^2 \alpha^{3/2}} [(D/D_0)^{1/3} - 1]^{1/2} \left[ \frac{D_0^{2/3}}{1 - D^{2/3}} \right] \qquad (21)$$

$$c_3 = \frac{2\sqrt{\alpha}}{\pi} [(D/D_0)^{1/3} - 1]^{1/2}. \qquad (22)$$

## 2.4. Weakening Function

We next turn to the formulation of the function $f(D)$ used to weaken the shear modulus as damage accumulates. The simplest effective medium representation of a cracked isotropic material assumes non-interacting cracks (Kachanov, 1993). Within this approximation, the change in elastic strain energy due to a population of cracks can be inferred by summing their individual contribution. This amounts to elastic compliances scaling linearly with damage. Elastic stiffnesses therefore scale as $(1 + CD)^{-1}$, where C is a constant that depends on the orientation distribution and geometry of cracks. Linearization of this form provides a reasonable estimate of elastic stiffnesses at low damage values. Because the damage framework of Ashby and Sammis (1990) sets an upper bound on damage at 1, we use this approximation and postulate a linear weakening of G with respect to damage D:

$$f(D) = \frac{\gamma - 1}{1 - D_0} D + \frac{1 - \gamma D_0}{1 - D_0}, \qquad (23)$$
the characteristic area of each wing crack ($\pi(l + \alpha a)^2$) amounts to the area $A_c$ that is obtained by projecting the spherical volume $V_c = 1/N_v$ along $\sigma_3$. Therefore, $A_b = A_c - \pi(l + \alpha a)^2$, with

$$A_c = \pi^{1/3} \left( \frac{3}{4N_v} \right)^{2/3}. \qquad (15)$$

This leads to the following expression for the internal stress acting in the direction of least compression, $\sigma_3^i$:

$$\sigma_3^i = \frac{F_w}{A_b}. \qquad (16)$$

Internal stress $\sigma_3^i$ increases dramatically as wings lengthen (D approaches 1) and the areas between fractures ($A_b$) shrink. This is when crack interactions become dominant. $K_I^{(i)}$ is readily obtained from $\sigma_3^i$ by analogy with Equation 14:

$$K_I^{(i)} = \frac{2}{\pi}(\sigma_3^i)\sqrt{\pi l}. \qquad (17)$$

The full expression of $K_I$ then reads

$$K_I = \frac{F_w}{[\pi(l + \beta a)]^{3/2}} + \frac{2}{\pi}(\sigma_3 + \sigma_3^i)\sqrt{\pi l}. \qquad (18)$$

It can be recast as a function of damage rather than crack length, following Bhat et al. (2011), yielding

$$K_I = \sqrt{\pi a}[(\sigma_3 A_3 - \sigma_1 A_1)(c_1 + c_2) + \sigma_3 c_3], \qquad (19)$$

where $c_1$, $c_2$, and $c_3$ are functions of the damage state that write

$$c_1 = \frac{1}{\pi^2 \alpha^{3/2} [(D/D_0)^{1/3} - 1 + \beta/\alpha]^{3/2}} \qquad (20)$$

$$c_2 = \frac{2}{\pi^2 \alpha^{3/2}} [(D/D_0)^{1/3} - 1]^{1/2} \left[ \frac{D_0^{2/3}}{1 - D^{2/3}} \right] \qquad (21)$$

$$c_3 = \frac{2\sqrt{\alpha}}{\pi} [(D/D_0)^{1/3} - 1]^{1/2}. \qquad (22)$$

## 2.4. Weakening Function

We next turn to the formulation of the function $f(D)$ used to weaken the shear modulus as damage accumulates. The simplest effective medium representation of a cracked isotropic material assumes non-interacting cracks (Kachanov, 1993). Within this approximation, the change in elastic strain energy due to a population of cracks can be inferred by summing their individual contribution. This amounts to elastic compliances scaling linearly with damage. Elastic stiffnesses therefore scale as $(1 + CD)^{-1}$, where C is a constant that depends on the orientation distribution and geometry of cracks. Linearization of this form provides a reasonable estimate of elastic stiffnesses at low damage values. Because the damage framework of Ashby and Sammis (1990) sets an upper bound on damage at 1, we use this approximation and postulate a linear weakening of G with respect to damage D:

$$f(D) = \frac{\gamma - 1}{1 - D_0} D + \frac{1 - \gamma D_0}{1 - D_0}, \qquad (23)$$
such that $G(D_0) = G_0$ and $G(D = 1) = \gamma G_0$. The weakening parameter $\gamma \in [0, 1]$ can be thought of as a property of the material representing the stiffness of a fully damaged rock (e.g., a fault zone) normalized by its maximum possible stiffness in a low-damage state. The derivative of our weakening function with respect to $D$ is:

$$f' = \frac{\gamma - 1}{1 - D_0} \qquad (24)$$

It should be noted that for simplicity our model weakens the shear modulus isotropically, even though damage grows in a highly anisotropic fashion.

To recap, Equations 4 and 5, combined with Equations 8, 19, and 23 make up the complete SCAM constitutive law, which is akin to Maxwell visco-elasticity with a strongly non linear dependence of viscosity on stress, and progressive alteration of the elastic modulus with increasing damage. These equations are reiterated below:

$$\begin{align*}
\dot{e}_{ij} &= \frac{\dot{s}_{ij}}{2G_0 f(D)} + \frac{s_{ij}}{2\eta_D} \\
\eta_D &= \frac{f^2(D) G_0}{|f'(D)| \dot{D}} \\
\dot{D} &= \frac{3D^2 D_0^{\frac{1}{4}} l_0}{\alpha a (K_I / K_{IC})^n} \\
K_I &= \sqrt{\pi a [(\sigma_3 A_3 - \sigma_1 A_1)(c_1 + c_2) + \sigma_3 c_3]} \\
f(D) &= \frac{\gamma - 1}{1 - D_0} D + \frac{1 - \gamma D_0}{1 - D_0}
\end{align*}$$

<figure><img src="image_6.png" /><figcaption>**Figure 4.** 0-D simulations based on the SCAM model in a typical triaxial setup with a confining pressure of 150 MPa (stress state shown in Panel a). (a and b) Differential stress and damage versus strain in a constant strain rate experiment ($10^{-5} s^{-1}$). (c and d) Deviatoric axial strain rate and damage versus time in an experiment where stress is kept constant after reaching the stress and damage state pictured by the stars in panels. (a) and (b) Here strains and stresses are represented positive in compression for clarity.</figcaption></figure>
such that $G(D_0) = G_0$ and $G(D = 1) = \gamma G_0$. The weakening parameter $\gamma \in [0, 1]$ can be thought of as a property of the material representing the stiffness of a fully damaged rock (e.g., a fault zone) normalized by its maximum possible stiffness in a low-damage state. The derivative of our weakening function with respect to $D$ is:

$$f' = \frac{\gamma - 1}{1 - D_0} \qquad (24)$$

It should be noted that for simplicity our model weakens the shear modulus isotropically, even though damage grows in a highly anisotropic fashion.

To recap, Equations 4 and 5, combined with Equations 8, 19, and 23 make up the complete SCAM constitutive law, which is akin to Maxwell visco-elasticity with a strongly non linear dependence of viscosity on stress, and progressive alteration of the elastic modulus with increasing damage. These equations are reiterated below:

$$\begin{align*}
\dot{e}_{ij} &= \frac{\dot{s}_{ij}}{2G_0 f(D)} + \frac{s_{ij}}{2\eta_D} \\
\eta_D &= \frac{f^2(D) G_0}{|f'(D)| \dot{D}} \\
\dot{D} &= \frac{3D^2 D_0^{\frac{1}{4}} l_0}{\alpha a (K_I / K_{IC})^n} \\
K_I &= \sqrt{\pi a [(\sigma_3 A_3 - \sigma_1 A_1)(c_1 + c_2) + \sigma_3 c_3]} \\
f(D) &= \frac{\gamma - 1}{1 - D_0} D + \frac{1 - \gamma D_0}{1 - D_0}
\end{align*}$$

<figure><img src="image_6.png" /><figcaption>**Figure 4.** 0-D simulations based on the SCAM model in a typical triaxial setup with a confining pressure of 150 MPa (stress state shown in Panel a). (a and b) Differential stress and damage versus strain in a constant strain rate experiment ($10^{-5} s^{-1}$). (c and d) Deviatoric axial strain rate and damage versus time in an experiment where stress is kept constant after reaching the stress and damage state pictured by the stars in panels. (a) and (b) Here strains and stresses are represented positive in compression for clarity.</figcaption></figure>
## 3. Application to a 0-D Triaxial Loading Setup

### 3.1. Constitutive SCAM Equations in a Triaxial Setup

To illustrate the behavior of the SCAM flow law, we implement it in a geometry typical of rock deformation experiments (Figure 4a): Compression along the axis of a cylindrical sample ($\sigma_1$, along direction $x_1$) subjected to axially symmetric confining stress ($\sigma_3$, along directions $x_2$ and $x_3$). The corresponding stress tensor writes

$$ \underline{\underline{\sigma}} = \begin{bmatrix} \sigma_{ax} & 0 & 0 \\ 0 & -p_c & 0 \\ 0 & 0 & -p_c \end{bmatrix}, \qquad (25) $$

where $\sigma_{ax}$ is the axial stress and $p_c$ the confining pressure surrounding the curved surface of the sample. We use a simplified point-wise formulation of our differential constitutive relationship Equation 4 assuming homogeneous deformation within the sample. As stated previously, we ignore volumetric strain and focus solely on the relationship between the deviatoric axial strain rate $\dot{e}_{ax}$ and the deviatoric axial stress $s_{ax}$. The constitutive equations reduce to the following ordinary differential equation (ODE):

$$ \dot{e}_{ax} = \frac{\dot{s}_{ax}}{2G_0 f(D)} + \frac{s_{ax}}{2\eta_D}, \qquad (26) $$

to be solved jointly with the damage evolution equation (Equation 8).

A first type of experiment consists of applying a constant axial strain rate and measuring the axial stress. In our framework, $s_{ax}$ verifies:

$$ \dot{s}_{ax} = 2G_0 f(D) \left( \dot{e}_{ax} - \frac{s_{ax}}{2\eta_D} \right), \qquad (27) $$

with $s_{ax}(t = 0) = 0$ and $D(t = 0) = D_i \ge D_0$.

Another class of experiments (brittle creep tests) consists of applying a constant axial stress and measuring the axial strain. In our model, the latter is given by

$$ \dot{e}_{ax} = \frac{s_{ax}}{2\eta_D}, \qquad (28) $$

In this case, the initial value of $D$ cannot be chosen arbitrarily and must be consistent with the imposed stress. To ensure this, we first integrate the constant strain rate and damage growth ODEs (Equations 8 and 27) up to the desired value of axial deviatoric stress $s_{ax}$ assuming a known strain rate. The damage value reached at the end of this preliminary step is used as initial condition for Equations 8 and 28, along with $e_{ax}(t = 0) = 0$. These equations are integrated up to $D$ close—but not equal—to 1, due to the singular behavior at this limit, coming from the $c_2$ term in Equation 21.

The above ODEs are integrated numerically using a 5/4th order Runge-Kutta method (Tsitouras et al., 2009). This is done within the DifferentialEquations.jl Julia package (Rackauckas & Nie, 2017) using adaptive time-stepping with absolute and relative tolerance of $10^{-6}$ and $10^{-4}$ respectively.

### 3.2. Stress-Strain Curves and Creep Regimes

We illustrate the fundamental behavior of the SCAM model in triaxial experiments using reference micro-mechanical parameters appropriate for Westerly granite (Table 1), which will be rigorously determined in Section 3.3. Constitutive equations are integrated up to $D = 0.95$.

Figures 4a and 4b correspond to a constant strain rate setup at $10^{-5}$ s$^{-1}$, under 150 MPa of confining pressure. The axial stress-strain curve displays an initial elastic phase followed by visible weakening of the effective modulus when differential stress exceeds ~700 MPa. This is accompanied by damage growth (Figure 4b) which
## 3. Application to a 0-D Triaxial Loading Setup

### 3.1. Constitutive SCAM Equations in a Triaxial Setup

To illustrate the behavior of the SCAM flow law, we implement it in a geometry typical of rock deformation experiments (Figure 4a): Compression along the axis of a cylindrical sample ($\sigma_1$, along direction $x_1$) subjected to axially symmetric confining stress ($\sigma_3$, along directions $x_2$ and $x_3$). The corresponding stress tensor writes

$$ \underline{\underline{\sigma}} = \begin{bmatrix} \sigma_{ax} & 0 & 0 \\ 0 & -p_c & 0 \\ 0 & 0 & -p_c \end{bmatrix}, \qquad (25) $$

where $\sigma_{ax}$ is the axial stress and $p_c$ the confining pressure surrounding the curved surface of the sample. We use a simplified point-wise formulation of our differential constitutive relationship Equation 4 assuming homogeneous deformation within the sample. As stated previously, we ignore volumetric strain and focus solely on the relationship between the deviatoric axial strain rate $\dot{e}_{ax}$ and the deviatoric axial stress $s_{ax}$. The constitutive equations reduce to the following ordinary differential equation (ODE):

$$ \dot{e}_{ax} = \frac{\dot{s}_{ax}}{2G_0 f(D)} + \frac{s_{ax}}{2\eta_D}, \qquad (26) $$

to be solved jointly with the damage evolution equation (Equation 8).

A first type of experiment consists of applying a constant axial strain rate and measuring the axial stress. In our framework, $s_{ax}$ verifies:

$$ \dot{s}_{ax} = 2G_0 f(D) \left( \dot{e}_{ax} - \frac{s_{ax}}{2\eta_D} \right), \qquad (27) $$

with $s_{ax}(t = 0) = 0$ and $D(t = 0) = D_i \ge D_0$.

Another class of experiments (brittle creep tests) consists of applying a constant axial stress and measuring the axial strain. In our model, the latter is given by

$$ \dot{e}_{ax} = \frac{s_{ax}}{2\eta_D}, \qquad (28) $$

In this case, the initial value of $D$ cannot be chosen arbitrarily and must be consistent with the imposed stress. To ensure this, we first integrate the constant strain rate and damage growth ODEs (Equations 8 and 27) up to the desired value of axial deviatoric stress $s_{ax}$ assuming a known strain rate. The damage value reached at the end of this preliminary step is used as initial condition for Equations 8 and 28, along with $e_{ax}(t = 0) = 0$. These equations are integrated up to $D$ close—but not equal—to 1, due to the singular behavior at this limit, coming from the $c_2$ term in Equation 21.

The above ODEs are integrated numerically using a 5/4th order Runge-Kutta method (Tsitouras et al., 2009). This is done within the DifferentialEquations.jl Julia package (Rackauckas & Nie, 2017) using adaptive time-stepping with absolute and relative tolerance of $10^{-6}$ and $10^{-4}$ respectively.

### 3.2. Stress-Strain Curves and Creep Regimes

We illustrate the fundamental behavior of the SCAM model in triaxial experiments using reference micro-mechanical parameters appropriate for Westerly granite (Table 1), which will be rigorously determined in Section 3.3. Constitutive equations are integrated up to $D = 0.95$.

Figures 4a and 4b correspond to a constant strain rate setup at $10^{-5}$ s$^{-1}$, under 150 MPa of confining pressure. The axial stress-strain curve displays an initial elastic phase followed by visible weakening of the effective modulus when differential stress exceeds ~700 MPa. This is accompanied by damage growth (Figure 4b) which
<table><caption>**Table 1**
Inverted Parameters</caption><thead><tr><th>Symbol</th><th>Description</th><th>Westerly granite</th><th>Darley Dale sandstone</th></tr></thead><tbody><tr><td>G<sub>0</sub></td><td>Shear modulus at D = D<sub>0</sub> (GPa)</td><td>28.72 ± 0.02</td><td>5.202 ± 3</td></tr><tr><td>γ</td><td>Residual ratio G/G<sub>0</sub> at D = 1</td><td>0.432 ± 0.004</td><td>0.281 ± 0.003</td></tr><tr><td>μ</td><td>Friction coefficient of the shear flaws</td><td>0.703 ± 0.001</td><td>0.5093 ± 0.0003</td></tr><tr><td>a</td><td>Shear flaws radius (μm)</td><td>6.66 ± 0.16</td><td>656.8 ± 4.3</td></tr><tr><td>n</td><td>Charles law exponent</td><td>11.82 ± 0.03</td><td>24.96 ± 0.07</td></tr><tr><td>i̇<sub>0</sub></td><td>Charles law reference crack growth rate (mm s<sup>-1</sup>)</td><td>16.44 ± 0.19</td><td>0.0029 ± 0.0001</td></tr><tr><td>K<sub>IC</sub></td><td>Fracture toughness (MPa m<sup>1/2</sup>)</td><td>1.29 ± 0.01</td><td>1.412 ± 0.004</td></tr><tr><td>D<sub>0</sub></td><td>D associated to the shear flaws only</td><td>0.1358 ± 0.0003</td><td>0.3724 ± 0.0008</td></tr><tr><td>D<sub>i</sub></td><td>Initial value of D</td><td>0.1361 ± 0.0003</td><td>0.27 ± 0.05</td></tr></tbody></table>

accelerates catastrophically as the sample reaches its peak stress. The post-peak stress drop is similarly abrupt as damage approaches 1.

Figures 4c and 4d correspond to a constant stress simulation starting at the yellow star shown in panels A and B. Strain rate first decelerates, then remains steady for hours, and ultimately accelerates up to the macroscopic failure of the material, consistent with the subsequent phases of brittle creep observed experimentally (Figure 1e). What is usually referred to as secondary creep would here be associated to the transition between decelerating and accelerating creep, and was not depicted in Figures 4c and 4d because of the clear bimodal dynamics of brittle creep expressed by the SCAM model. This strain rate behavior is associated with dynamics similar to those of the damage growth rate, visible through the slope of damage evolution with respect to time in Figure 4d.

The effect of various model parameters and experimental conditions on the behavior of the SCAM model under constant strain rate is shown in Figure 5. The black curves correspond to a strain rate of $10^{-5}$ s$^{-1}$, a confining pressure of 150 MPa and the reference set of micromechanical parameters for Westerly granite (Table 1). Figures 5a and 5b shows that a reduction of the imposed axial strain rate leads to a lower peak stress due to damage having more time to accumulate under lower axial stress, precipitating failure (Figure 5b). Increasing the radius of the shear defects (Panel e) while keeping $D_0$ constant leads to a decrease of the peak stress. This is because the stress intensity factor at the wing crack tips increases with increasing shear defect size, prompting faster crack growth. Thus, significant damage can build under lower stresses, and the peak stress is reached sooner. Increasing $D_0$ while keeping the shear defect size constant (Panel f) also leads to a lower peak stress, but limits the amount of softening that takes place pre-peak. This is because cracks arranged in a denser array will interact and coalesce sooner. The stress decrease additionally does not display the abrupt drop seen with the reference case, which we attribute to the larger reduction of shear modulus per damage increment. A decrease of the Charles law exponent $n$ (Panel g) or friction coefficient (Panel d) similarly lowers the peak stress, by enabling damage build-up under lower stress intensity factors Equation 7 and under lower differential stress, respectively. Finally, a greater degree of modulus weakening (via a reduced $γ$ parameter) leads to more pre-peak softening but has a limited impact on the peak stress (Panel h). It also stabilizes the stress drop by limiting the unstable growth of damage as it gets close to 1. In this case of extreme loss of elastic stiffness, the larger negative stress increment associated with damage increment post-peak tends to reduce the catastrophic increase in damage growth.

To better visualize the dynamics of damage growth in the SCAM model, we represent constant strain rate experiments in a plot of differential stress versus damage (Figure 6a). This representation allows us to map the stress intensity factor at the wing-crack tip (colors and contours in panel a), which gives us a proxy for damage growth rate. We specifically highlight two sets of experiments. The first set is performed at a laboratory strain rate $\dot{e}_{ax} = 10^{-5}$ s$^{-1}$, and the second at a tectonic strain rate of $10^{-15}$ s$^{-1}$, both under a confining pressure $p_c = 150$ MPa. In each set, we vary the initial damage $D_i$, using values of 0.136 ($D_0$), 0.227, 0.318, 0.409, and 0.5.

Each experiment follows a specific trajectory in differential stress versus damage space. For example, in the case of no initial tensile cracks ($D_i = D_0$), differential stress first increases while damage remains constant. This is because in the initial elastic regime, $K_I \le 0$ and wing cracks cannot grow. Once the system reaches the domain of positive $K_I$, damage can start growing, and increases with stress. The system appears to follow a contour of
<table><caption>**Table 1**
Inverted Parameters</caption><thead><tr><th>Symbol</th><th>Description</th><th>Westerly granite</th><th>Darley Dale sandstone</th></tr></thead><tbody><tr><td>G<sub>0</sub></td><td>Shear modulus at D = D<sub>0</sub> (GPa)</td><td>28.72 ± 0.02</td><td>5.202 ± 3</td></tr><tr><td>γ</td><td>Residual ratio G/G<sub>0</sub> at D = 1</td><td>0.432 ± 0.004</td><td>0.281 ± 0.003</td></tr><tr><td>μ</td><td>Friction coefficient of the shear flaws</td><td>0.703 ± 0.001</td><td>0.5093 ± 0.0003</td></tr><tr><td>a</td><td>Shear flaws radius (μm)</td><td>6.66 ± 0.16</td><td>656.8 ± 4.3</td></tr><tr><td>n</td><td>Charles law exponent</td><td>11.82 ± 0.03</td><td>24.96 ± 0.07</td></tr><tr><td>i̇<sub>0</sub></td><td>Charles law reference crack growth rate (mm s<sup>-1</sup>)</td><td>16.44 ± 0.19</td><td>0.0029 ± 0.0001</td></tr><tr><td>K<sub>tc</sub></td><td>Fracture toughness (MPa m<sup>1/2</sup>)</td><td>1.29 ± 0.01</td><td>1.412 ± 0.004</td></tr><tr><td>D<sub>0</sub></td><td>D associated to the shear flaws only</td><td>0.1358 ± 0.0003</td><td>0.3724 ± 0.0008</td></tr><tr><td>D<sub>i</sub></td><td>Initial value of D</td><td>0.1361 ± 0.0003</td><td>0.27 ± 0.05</td></tr></tbody></table>

accelerates catastrophically as the sample reaches its peak stress. The post-peak stress drop is similarly abrupt as damage approaches 1.

Figures 4c and 4d correspond to a constant stress simulation starting at the yellow star shown in panels A and B. Strain rate first decelerates, then remains steady for hours, and ultimately accelerates up to the macroscopic failure of the material, consistent with the subsequent phases of brittle creep observed experimentally (Figure 1e). What is usually referred to as secondary creep would here be associated to the transition between decelerating and accelerating creep, and was not depicted in Figures 4c and 4d because of the clear bimodal dynamics of brittle creep expressed by the SCAM model. This strain rate behavior is associated with dynamics similar to those of the damage growth rate, visible through the slope of damage evolution with respect to time in Figure 4d.

The effect of various model parameters and experimental conditions on the behavior of the SCAM model under constant strain rate is shown in Figure 5. The black curves correspond to a strain rate of $10^{-5}$ s$^{-1}$, a confining pressure of 150 MPa and the reference set of micromechanical parameters for Westerly granite (Table 1). Figures 5a and 5b shows that a reduction of the imposed axial strain rate leads to a lower peak stress due to damage having more time to accumulate under lower axial stress, precipitating failure (Figure 5b). Increasing the radius of the shear defects (Panel e) while keeping $D_0$ constant leads to a decrease of the peak stress. This is because the stress intensity factor at the wing crack tips increases with increasing shear defect size, prompting faster crack growth. Thus, significant damage can build under lower stresses, and the peak stress is reached sooner. Increasing $D_0$ while keeping the shear defect size constant (Panel f) also leads to a lower peak stress, but limits the amount of softening that takes place pre-peak. This is because cracks arranged in a denser array will interact and coalesce sooner. The stress decrease additionally does not display the abrupt drop seen with the reference case, which we attribute to the larger reduction of shear modulus per damage increment. A decrease of the Charles law exponent $n$ (Panel g) or friction coefficient (Panel d) similarly lowers the peak stress, by enabling damage build-up under lower stress intensity factors Equation 7 and under lower differential stress, respectively. Finally, a greater degree of modulus weakening (via a reduced $γ$ parameter) leads to more pre-peak softening but has a limited impact on the peak stress (Panel h). It also stabilizes the stress drop by limiting the unstable growth of damage as it gets close to 1. In this case of extreme loss of elastic stiffness, the larger negative stress increment associated with damage increment post-peak tends to reduce the catastrophic increase in damage growth.

To better visualize the dynamics of damage growth in the SCAM model, we represent constant strain rate experiments in a plot of differential stress versus damage (Figure 6a). This representation allows us to map the stress intensity factor at the wing-crack tip (colors and contours in panel a), which gives us a proxy for damage growth rate. We specifically highlight two sets of experiments. The first set is performed at a laboratory strain rate $\dot{e}_{ax} = 10^{-5}$ s$^{-1}$, and the second at a tectonic strain rate of $10^{-15}$ s$^{-1}$, both under a confining pressure $p_c = 150$ MPa. In each set, we vary the initial damage $D_i$, using values of 0.136 ($D_0$), 0.227, 0.318, 0.409, and 0.5.

Each experiment follows a specific trajectory in differential stress versus damage space. For example, in the case of no initial tensile cracks ($D_i = D_0$), differential stress first increases while damage remains constant. This is because in the initial elastic regime, $K_I \le 0$ and wing cracks cannot grow. Once the system reaches the domain of positive $K_I$, damage can start growing, and increases with stress. The system appears to follow a contour of
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Effect of strain rate (a), (b), confining pressure $p_c$ (c), friction coefficient $\mu$ (d), shear defect radius $a$ (e), shear crack density $D_0$ (f), corrosion index $n$ (g) and residual shear modulus $\gamma$ (h) on differential stress with respect to deviatoric axial strain for 0-D SCAM simulations. Black lines correspond to the best fitting parameters for Westerly granite detailed in Section 3.3. Here strains and stresses are represented positive in compression for clarity.</figcaption></figure>

constant $K_I$ up to the peak differential stress (~1,140 MPa). Past this point, the differential stress starts to decrease while damage keeps increasing at an accelerating pace. This is due to the fact that $K_I$, which sets the rate of damage growth, now increases with increasing damage. This final phase of rapid failure manifests as an abrupt post-peak stress drop in the stress-strain curve (Panel b).

These three regimes, characterized by the absence of growth, the stable and then the unstable growth of damage is illustrated in Figure 3 with the three numbered stars respectively. Simulations carried out under the same strain rate, but with greater initial damage show the same behavior, and their trajectories tend to align along the same iso-$K_I$ (~0.3$K_{IC}$) path as followed by the $D_I = D_0$ case. This forms an envelope that materializes an upper bound of the differential stress value with respect to damage. This envelope corresponds to ($K_I$ ~ 0.05$K_{IC}$) for tectonic strain rates, and therefore lies at lower stress values. If, however, a simulation is initiated with damage in excess of ~0.45 (e.g., orange paths in panel A), damage will immediately start growing in the unstable regime, where $\partial K_I / \partial D > 0$. In this case, the system reaches a peak stress which is lower than that of the other simulations.
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Effect of strain rate (a), (b), confining pressure $p_c$ (c), friction coefficient $\mu$ (d), shear defect radius $a$ (e), shear crack density $D_0$ (f), corrosion index $n$ (g) and residual shear modulus $\gamma$ (h) on differential stress with respect to deviatoric axial strain for 0-D SCAM simulations. Black lines correspond to the best fitting parameters for Westerly granite detailed in Section 3.3. Here strains and stresses are represented positive in compression for clarity.</figcaption></figure>

constant $K_I$ up to the peak differential stress (~1,140 MPa). Past this point, the differential stress starts to decrease while damage keeps increasing at an accelerating pace. This is due to the fact that $K_I$, which sets the rate of damage growth, now increases with increasing damage. This final phase of rapid failure manifests as an abrupt post-peak stress drop in the stress-strain curve (Panel b).

These three regimes, characterized by the absence of growth, the stable and then the unstable growth of damage is illustrated in Figure 3 with the three numbered stars respectively. Simulations carried out under the same strain rate, but with greater initial damage show the same behavior, and their trajectories tend to align along the same iso-$K_I$ (~0.3$K_{IC}$) path as followed by the $D_I = D_0$ case. This forms an envelope that materializes an upper bound of the differential stress value with respect to damage. This envelope corresponds to ($K_I$ ~ 0.05$K_{IC}$) for tectonic strain rates, and therefore lies at lower stress values. If, however, a simulation is initiated with damage in excess of ~0.45 (e.g., orange paths in panel A), damage will immediately start growing in the unstable regime, where $\partial K_I / \partial D > 0$. In this case, the system reaches a peak stress which is lower than that of the other simulations.
<figure><img src="image_8.png" /><figcaption>**Figure 6.** (a) Trajectories of differential stress with respect to damage at $\dot{e}_{ax} = -10^{-5} s^{-1}$ and $\dot{e}_{ax} = -10^{-15} s^{-1}$ for various initial damages (color code) and at constant confining pressure $\rho_c = 150$ MPa. The background is colored according to $K_\phi/K_{IC}$, with the two iso-values pictured with black dashed contours representing $K_I = 0$ and $K_I = K_{IC}$ respectively. The near-vertical dashed black line highlights damage values where $\partial K_\phi/\partial D = 0$ under constant $\Delta\sigma$. (b) Same trajectories plotted as standard deviatoric axial strain versus differential stress, along with the color code for initial damage state. Strains and stresses are represented positive in compression for clarity.</figcaption></figure>

### 3.3. Calibration of SCAM Parameters With Laboratory Experiments

The ability of the SCAM model to reproduce both constant strain rate and constant stress experiments suggests that laboratory data can be used to constrain its micromechanical parameters (Table 1). Specifically, stress-strain curves from constant strain rate experiments under various confining pressures can help constrain elastic and frictional properties, while strain rates and time to failure in brittle creep tests contain information about the kinetics of damage build-up.

To leverage this information, we use the 0-D “forward” models presented in the previous section in a Bayesian inversion framework (Tarantola, 2005, see Appendix A for details). We expect 0-D models to be representative of the homogeneous deformation stage up to the peak stress (prior to localization), as micro-cracking is known to first develop in a distributed fashion (Figure 1).

#### 3.3.1. Experimental Data

We apply the Bayesian inversion method to experimental data corresponding to two lithologies. The first is Westerly granite, a rock type widely used in experiments that is representative of the continental upper crust in term of mineralogy and low porosity. This rock has been shown to experience the type of diffuse cracking and catastrophic fracture coalescence that our model seeks to capture (e.g., Lockner et al., 1991; Tapponnier & Brace, 1976). We specifically use constant strain rate ($\dot{e}_{ax} = 10^{-5} s^{-1}$) experiments under confining pressures of 20, 30, 80, and 150 MPa in dry conditions from Wawersik and Brace (1971) (Figure 7a). We complement these data with minimum brittle (secondary) creep strain rates measured under seven imposed differential stresses ranging from 77% to 93% of the short-term strength (meaning the peak strength at a laboratory strain rate) of the
<figure><img src="image_8.png" /><figcaption>**Figure 6.** (a) Trajectories of differential stress with respect to damage at $\dot{e}_{ax} = -10^{-5} s^{-1}$ and $\dot{e}_{ax} = -10^{-15} s^{-1}$ for various initial damages (color code) and at constant confining pressure $\rho_c = 150$ MPa. The background is colored according to $K_\phi/K_{IC}$, with the two iso-values pictured with black dashed contours representing $K_I = 0$ and $K_I = K_{IC}$ respectively. The near-vertical dashed black line highlights damage values where $\partial K_\phi/\partial D = 0$ under constant $\Delta\sigma$. (b) Same trajectories plotted as standard deviatoric axial strain versus differential stress, along with the color code for initial damage state. Strains and stresses are represented positive in compression for clarity.</figcaption></figure>

### 3.3. Calibration of SCAM Parameters With Laboratory Experiments

The ability of the SCAM model to reproduce both constant strain rate and constant stress experiments suggests that laboratory data can be used to constrain its micromechanical parameters (Table 1). Specifically, stress-strain curves from constant strain rate experiments under various confining pressures can help constrain elastic and frictional properties, while strain rates and time to failure in brittle creep tests contain information about the kinetics of damage build-up.

To leverage this information, we use the 0-D “forward” models presented in the previous section in a Bayesian inversion framework (Tarantola, 2005, see Appendix A for details). We expect 0-D models to be representative of the homogeneous deformation stage up to the peak stress (prior to localization), as micro-cracking is known to first develop in a distributed fashion (Figure 1).

#### 3.3.1. Experimental Data

We apply the Bayesian inversion method to experimental data corresponding to two lithologies. The first is Westerly granite, a rock type widely used in experiments that is representative of the continental upper crust in term of mineralogy and low porosity. This rock has been shown to experience the type of diffuse cracking and catastrophic fracture coalescence that our model seeks to capture (e.g., Lockner et al., 1991; Tapponnier & Brace, 1976). We specifically use constant strain rate ($\dot{e}_{ax} = 10^{-5} s^{-1}$) experiments under confining pressures of 20, 30, 80, and 150 MPa in dry conditions from Wawersik and Brace (1971) (Figure 7a). We complement these data with minimum brittle (secondary) creep strain rates measured under seven imposed differential stresses ranging from 77% to 93% of the short-term strength (meaning the peak strength at a laboratory strain rate) of the
<figure><img src="image_9.png" /><figcaption>**Figure 7.** SCAM model fit to experimental data on Westerly granite. (a) Constant strain rate experiments with $\dot{e}_{ax} = -10^{-5} \text{ s}^{-1}$ and confining pressures of 20 (orange), 30 (red), 80 (blue) and 150 MPa (black). Dots show the data and lines the best fitting models. (b) Dataset of minimum brittle creep strain rate as a function of imposed differential stress (black dots), with relationship derived from best fitting model (red line). (c) Prior (gray) and posterior (black) distributions of inverted micromechanical parameters.</figcaption></figure>

rock subjected to an effective confining pressure of 30 MPa in water-saturated samples by Brantut et al. (2012) (Figure 7b).

The second rock type we consider is Darley Dale feldspar-rich sandstone, another widely studied lithology. While its properties are likely less representative of the upper crust than that of Westerly granite, inferring its micro-mechanical parameters can provide helpful comparisons to assess the validity of our model. One caveat of this choice is that porous sandstone may deform according to mechanisms other than the growth of tensile cracks from shear defects, such as Hertz-contact driven microfracturing (e.g., Zhang et al., 1990), tensile cracks nucleating from pores (e.g., Sammis & Ashby, 1986), or distributed cataclastic flow where microcracking grows along grain boundaries (e.g., Menéndez et al., 1996). Heap et al. (2009) however report that for confining pressures up to their maximum of 50 MPa, stress-induced damage grows predominantly in a direction subparallel to the axis of compression. Additionally, dilatancy patterns at constant strain rate are very similar to what is observed in low-porosity rocks such as Westerly granite (e.g., X. Wu et al., 2000; Zoback & Byerlee, 1975). These observations suggest that the wing-crack model remains a relevant conceptual framework for the brittle deformation of Darley Dale sandstone, at least below 50 MPa of effective confining pressure.
<figure><img src="image_9.png" /><figcaption>**Figure 7.** SCAM model fit to experimental data on Westerly granite. (a) Constant strain rate experiments with $\dot{e}_{ax} = -10^{-5} \text{ s}^{-1}$ and confining pressures of 20 (orange), 30 (red), 80 (blue) and 150 MPa (black). Dots show the data and lines the best fitting models. (b) Dataset of minimum brittle creep strain rate as a function of imposed differential stress (black dots), with relationship derived from best fitting model (red line). (c) Prior (gray) and posterior (black) distributions of inverted micromechanical parameters.</figcaption></figure>

rock subjected to an effective confining pressure of 30 MPa in water-saturated samples by Brantut et al. (2012) (Figure 7b).

The second rock type we consider is Darley Dale feldspar-rich sandstone, another widely studied lithology. While its properties are likely less representative of the upper crust than that of Westerly granite, inferring its micro-mechanical parameters can provide helpful comparisons to assess the validity of our model. One caveat of this choice is that porous sandstone may deform according to mechanisms other than the growth of tensile cracks from shear defects, such as Hertz-contact driven microfracturing (e.g., Zhang et al., 1990), tensile cracks nucleating from pores (e.g., Sammis & Ashby, 1986), or distributed cataclastic flow where microcracking grows along grain boundaries (e.g., Menéndez et al., 1996). Heap et al. (2009) however report that for confining pressures up to their maximum of 50 MPa, stress-induced damage grows predominantly in a direction subparallel to the axis of compression. Additionally, dilatancy patterns at constant strain rate are very similar to what is observed in low-porosity rocks such as Westerly granite (e.g., X. Wu et al., 2000; Zoback & Byerlee, 1975). These observations suggest that the wing-crack model remains a relevant conceptual framework for the brittle deformation of Darley Dale sandstone, at least below 50 MPa of effective confining pressure.
Friction coefficient $\mu$, fracture toughness $K_{IC}$ and pre-existing crack length $a$ were previously estimated for Darley-Dale by X. Wu et al. (2000) within the Ashby and Sammis (1990) wing-crack framework, but assuming critical crack growth (i.e., at $K_I = K_{IC}$). Here, we additionally use entire time series of brittle creep strain rates from Heap et al. (2009) to provide strong constraints on the kinetics of crack growth. We specifically use experimental results performed under constant stresses of 80%, 85% and 90% of the short-term strength from Heap et al. (2009) (Figure 8b). We combine these time series with stress-strain curves determined under a constant strain rate of $10^{-5} s^{-1}$ and confining pressures of 10 and 50 MPa in water-saturated samples from the same authors (Figure 8a). While data was also available for a confining pressure of 30 MPa, it displayed a significantly different shear modulus compared to the other two. We thus decided to exclude it from the joint inversion procedure.

It should be noted that because our constitutive law is based on incompressible elasticity, we remove the volumetric component of elastic deformation (i.e., the Poisson effect) from constant strain rate experimental data. In

<figure><img src="image_10.png" /><figcaption>**Figure 8.** SCAM model fit to experimental data on Darley Dale sandstone. (a) Constant strain rate experiments with $\dot{e}_{ax} = -10^{-5} s^{-1}$ and confining pressures of 10 (orange) and 50 MPa (blue). Dots show the data and lines the best fitting models. Larger circles and squares indicate experimentally determined and modeled peak stress, respectively. Horizontal dashed lines show the imposed stress levels for the experiments shown in Panel (b). (b) Brittle creep tests shown as time series of axial strain rate under imposed differential stress and confining pressure $p_c = 30$ MPa (dots and lines correspond to data and models, respectively). A simulated constant strain rate curve at this confining pressure is included in Panel a (red dashed curve with red square marking the peak stress). (c) Prior (gray) and posterior (black) distributions of inverted parameters. The posterior distribution of the initial damage state $D_i$ is the same as its prior because the inversion kept lowering the $D_i$ value to below $D_0$, which was not permitted. This resulted in a null gradient of the log-likelihood with respect to $D_i$ and thus no change of the parameter nor of its posterior distribution relative to the prior.</figcaption></figure>
Friction coefficient $\mu$, fracture toughness $K_{IC}$ and pre-existing crack length $a$ were previously estimated for Darley-Dale by X. Wu et al. (2000) within the Ashby and Sammis (1990) wing-crack framework, but assuming critical crack growth (i.e., at $K_I = K_{IC}$). Here, we additionally use entire time series of brittle creep strain rates from Heap et al. (2009) to provide strong constraints on the kinetics of crack growth. We specifically use experimental results performed under constant stresses of 80%, 85% and 90% of the short-term strength from Heap et al. (2009) (Figure 8b). We combine these time series with stress-strain curves determined under a constant strain rate of $10^{-5} s^{-1}$ and confining pressures of 10 and 50 MPa in water-saturated samples from the same authors (Figure 8a). While data was also available for a confining pressure of 30 MPa, it displayed a significantly different shear modulus compared to the other two. We thus decided to exclude it from the joint inversion procedure.

It should be noted that because our constitutive law is based on incompressible elasticity, we remove the volumetric component of elastic deformation (i.e., the Poisson effect) from constant strain rate experimental data. In

<figure><img src="image_10.png" /><figcaption>**Figure 8.** SCAM model fit to experimental data on Darley Dale sandstone. (a) Constant strain rate experiments with $\dot{e}_{ax} = -10^{-5} s^{-1}$ and confining pressures of 10 (orange) and 50 MPa (blue). Dots show the data and lines the best fitting models. Larger circles and squares indicate experimentally determined and modeled peak stress, respectively. Horizontal dashed lines show the imposed stress levels for the experiments shown in Panel (b). (b) Brittle creep tests shown as time series of axial strain rate under imposed differential stress and confining pressure $p_c = 30$ MPa (dots and lines correspond to data and models, respectively). A simulated constant strain rate curve at this confining pressure is included in Panel a (red dashed curve with red square marking the peak stress). (c) Prior (gray) and posterior (black) distributions of inverted parameters. The posterior distribution of the initial damage state $D_i$ is the same as its prior because the inversion kept lowering the $D_i$ value to below $D_0$, which was not permitted. This resulted in a null gradient of the log-likelihood with respect to $D_i$ and thus no change of the parameter nor of its posterior distribution relative to the prior.</figcaption></figure>
practice, this means that the reported strain rate is converted to a deviatoric axial strain rate $\dot{\epsilon}_{ax}' = \frac{2}{3}(1 + \nu)\dot{\epsilon}_{ax}$, with $\nu$ the Poisson's ratio, in order to compare our simulations with laboratory results.

### 3.3.2. A-Priori Parameters

The a-priori value of $G$ is chosen based on the initial slope of the elastic (linear) portion of available stress-strain curves. The initial value of the friction coefficient $\mu$ is set to the standard value of 0.65 (Byerlee, 1978). We initially assume $K_{IC}$ to fall between values determined experimentally in quartz (1 MPa $\cdot m^{1/2}$) and wet Westerly granite (1.74 MPa $\cdot m^{1/2}$) from the compilation of Atkinson (1984). We initialize $\gamma$ at an intermediate value of 0.5, and $a$ at the mean grain size of Westerly granite: 0.5 mm. We also set $D_0 = D_i$ at 0.2, $n = 12$ (Atkinson, 1979), and $\dot{l}_0 = 10^{-2} \text{ m s}^{-1}$. This set of a-priori guesses on the parameter values is first used to invert only the data from constant strain rate experiments. The results of this step are used to construct new priors on the model parameters, shown as gray shadings in the right columns of Figures 7 and 8. These priors are then used for a combined inversion of constant strain rate and brittle creep experiments, using a step multiplier $\kappa = 0.1$ (See Appendix A for details). A hundred steps were typically sufficient to reach convergence, yielding the posterior model parameter distributions shown in black in the right column of Figures 7 and 8.

### 3.3.3. Results

Figure 7 shows the results of our joint inversion of constant strain rate and brittle creep data in Westerly granite. Panel a compares the SCAM-simulated stress-strain curves (plain lines) and the experimental data points. The agreement is good at confining pressures of 30, 80, and 150 MPa. At 20 MPa, however, the model slightly overestimates the peak stress. The Wawersik and Brace (1971) study also contains data at atmospheric, 3.5 and 10 MPa confining pressures, but under these conditions our model was not able to accurately represent the pressure dependence of the peak stress. Figure 7b also shows secondary creep strain rates as a function of imposed differential stress from our simulations (red line), which are in good agreement with experimental data (black dots). The relationship between brittle creep strain rate and stress is effectively a power law with a stress exponent of ~18.

The best fitting parameter values as well as their log-normal standard deviations are listed in Table 1, and shown as probability distributions (in black) in Figure 7. Because the prior distributions (in gray) were determined by fitting only constant strain rate data, the differences between prior and posterior distributions highlight the information provided by brittle creep data. This information specifically constrains the initial damage state, as well as parameters related to the kinetics of damage growth such as $\dot{l}_0$ or $n$. It also strongly constrains the size of shear defects (to ~7 µm), which influences $K_I$ and therefore the damage growth rate. Parameters such as $K_{IC}$, $\mu$ and $\gamma$ are also slightly re-evaluated.

Best-fitting stress-strain curves for Darley Dale sandstone are shown in Figure 8a for 10 and 50 MPa of confining pressure as plain lines, along with experimental data (points). Larger markers mark the peak stress of simulations (squares) and experiments (circles). The red dashed line shows an additional simulated curve at an intermediate pressure of 30 MPa. This pressure corresponds to that of the brittle creep tests (Figure 8b), which were conducted under three axial stresses indicated as dashed lines in Figure 8a. It can be seen that the greatest applied differential stress (141 MPa) is very close to the inferred peak stress at 30 MPa of confining pressure. Figure 8b compares simulated and measured strain rates in the brittle creep experiments. Our best fitting parameters do a good job at reproducing the shape of the strain rate curves as well as the time to macroscopic failure (the final, near-vertical portion of the curves).

Similarly to our results in Westerly granite, joint inversions of brittle creep tests and constant strain rate experiments provide strong constraints on parameters such as shear defect size, initial damage, Charles law exponent, and $\dot{l}_0$. Our inversions yield a significantly greater defect size (~700 µm vs. ~7 µm) and Charles law exponent (~25 vs. ~12) in sandstone compared to granite, as well as a lower shear modulus and greater degree of elastic weakening (lower $\gamma$). The initial damage state of sandstone also appears greater. We however find comparable fracture toughness in both lithologies, and a slightly greater coefficient of (defect-scale) friction in granite (0.7 vs. 0.5).
practice, this means that the reported strain rate is converted to a deviatoric axial strain rate $\dot{\epsilon}_{ax}' = \frac{2}{3}(1 + \nu)\dot{\epsilon}_{ax}$, with $\nu$ the Poisson's ratio, in order to compare our simulations with laboratory results.

### 3.3.2. A-Priori Parameters

The a-priori value of $G$ is chosen based on the initial slope of the elastic (linear) portion of available stress-strain curves. The initial value of the friction coefficient $\mu$ is set to the standard value of 0.65 (Byerlee, 1978). We initially assume $K_{IC}$ to fall between values determined experimentally in quartz (1 MPa $\cdot m^{1/2}$) and wet Westerly granite (1.74 MPa $\cdot m^{1/2}$) from the compilation of Atkinson (1984). We initialize $\gamma$ at an intermediate value of 0.5, and $a$ at the mean grain size of Westerly granite: 0.5 mm. We also set $D_0 = D_i$ at 0.2, $n = 12$ (Atkinson, 1979), and $\dot{l}_0 = 10^{-2} \text{ m s}^{-1}$. This set of a-priori guesses on the parameter values is first used to invert only the data from constant strain rate experiments. The results of this step are used to construct new priors on the model parameters, shown as gray shadings in the right columns of Figures 7 and 8. These priors are then used for a combined inversion of constant strain rate and brittle creep experiments, using a step multiplier $\kappa = 0.1$ (See Appendix A for details). A hundred steps were typically sufficient to reach convergence, yielding the posterior model parameter distributions shown in black in the right column of Figures 7 and 8.

### 3.3.3. Results

Figure 7 shows the results of our joint inversion of constant strain rate and brittle creep data in Westerly granite. Panel a compares the SCAM-simulated stress-strain curves (plain lines) and the experimental data points. The agreement is good at confining pressures of 30, 80, and 150 MPa. At 20 MPa, however, the model slightly overestimates the peak stress. The Wawersik and Brace (1971) study also contains data at atmospheric, 3.5 and 10 MPa confining pressures, but under these conditions our model was not able to accurately represent the pressure dependence of the peak stress. Figure 7b also shows secondary creep strain rates as a function of imposed differential stress from our simulations (red line), which are in good agreement with experimental data (black dots). The relationship between brittle creep strain rate and stress is effectively a power law with a stress exponent of ~18.

The best fitting parameter values as well as their log-normal standard deviations are listed in Table 1, and shown as probability distributions (in black) in Figure 7. Because the prior distributions (in gray) were determined by fitting only constant strain rate data, the differences between prior and posterior distributions highlight the information provided by brittle creep data. This information specifically constrains the initial damage state, as well as parameters related to the kinetics of damage growth such as $\dot{l}_0$ or $n$. It also strongly constrains the size of shear defects (to ~7 µm), which influences $K_I$ and therefore the damage growth rate. Parameters such as $K_{IC}$, $\mu$ and $\gamma$ are also slightly re-evaluated.

Best-fitting stress-strain curves for Darley Dale sandstone are shown in Figure 8a for 10 and 50 MPa of confining pressure as plain lines, along with experimental data (points). Larger markers mark the peak stress of simulations (squares) and experiments (circles). The red dashed line shows an additional simulated curve at an intermediate pressure of 30 MPa. This pressure corresponds to that of the brittle creep tests (Figure 8b), which were conducted under three axial stresses indicated as dashed lines in Figure 8a. It can be seen that the greatest applied differential stress (141 MPa) is very close to the inferred peak stress at 30 MPa of confining pressure. Figure 8b compares simulated and measured strain rates in the brittle creep experiments. Our best fitting parameters do a good job at reproducing the shape of the strain rate curves as well as the time to macroscopic failure (the final, near-vertical portion of the curves).

Similarly to our results in Westerly granite, joint inversions of brittle creep tests and constant strain rate experiments provide strong constraints on parameters such as shear defect size, initial damage, Charles law exponent, and $\dot{l}_0$. Our inversions yield a significantly greater defect size (~700 µm vs. ~7 µm) and Charles law exponent (~25 vs. ~12) in sandstone compared to granite, as well as a lower shear modulus and greater degree of elastic weakening (lower $\gamma$). The initial damage state of sandstone also appears greater. We however find comparable fracture toughness in both lithologies, and a slightly greater coefficient of (defect-scale) friction in granite (0.7 vs. 0.5).
# 4. Application to a 2-D Plane-Strain Numerical Press

## 4.1. Conservation Equations and Numerical Methodology

In order to perform 2-D simulations of material deformation governed by the SCAM constitutive equations, we adapt the long-term tectonic modeling code SiStER (Simple Stokes solver with Exotic Rheologies, Olive et al. (2016)), which solves for conservation of mass, momentum (and energy if needed), in a 2-D continuum assuming elastic incompressibility and planar deformation. Conservation of mass and momentum write:

$$ \frac{\partial v_i}{\partial x_i} = 0, \qquad (29) $$

and

$$ \frac{\partial s_{ij}}{\partial x_j} - \frac{\partial P}{\partial x_i} + \rho g_i = 0, \qquad (30) $$

where $v_i$ are velocities, $P = -1/3 \sigma_{kk}$ is pressure, $\rho$ is density and $g_i$ the gravitational acceleration. Deviatoric stresses $s_{ij}$ are related to velocities in Equation 30 using a Maxwell visco-elastic constitutive relationship between deviatoric stresses and strain rates (e.g., Gerya, 2010; Moresi et al., 2003):

$$ \dot{e}_{ij} = \frac{1}{2G} \dot{s}_{ij} + \frac{1}{2\eta} s_{ij} \qquad (31) $$

where $\dot{s}_{ij}$ is discretized using a first-order backward finite difference scheme with time step $\Delta t$, so that the deviatoric stress at time $t$ becomes

$$ s_{ij}^t = 2Z\eta\dot{e}_{ij}^t + (1-Z)s_{ij}^{t-\Delta t}, \qquad (32) $$

with $\dot{e}_{ij} = (\partial v_i / \partial x_j + \partial v_j / \partial x_i) / 2$ and

$$ Z = \frac{G\Delta t}{\eta + G\Delta t}. \qquad (33) $$

The effective viscosity in Equation 31 can represent a range of rheologies. A very high value sets a very long Maxwell time, which effectively renders the material elastic. In the viscous regime, $\eta$ can represent brittle plasticity (e.g., as detailed in Section 4.2.2), or a specific creep mechanism of known flow law. In practice, $\eta$ is constructed as the harmonic average of several viscosities, each representing individual flow mechanisms.

The mass Equation 29 and momentum Equation 30 conservation equations, expressed in terms of velocities Equation 32 and pressure, are discretized with a conservative finite difference scheme formulated on a staggered grid (e.g., Gerya & Yuen, 2003). This leads to a linear system that is solved for velocities and pressure over the entire domain using a direct solver. Retroactions between the viscosity and velocity fields require the use of non-linear iterations (here approximate-Newton, described as Algorithm 2 in Spiegelman et al. (2016)) to reach convergence, which is assessed by comparing the L2 norm of the residual vector to a specified tolerance (relative tolerance between $10^{-7}$ and $10^{-2}$, see readme documentation in the code repository linked in the Acknowledgments section).

Once a reasonably converged solution is found, the time evolution is performed explicitly using the time step $\Delta t$ introduced in Equation 32. This is specifically done by advecting Lagrangian markers which carry material properties such as density, viscosity and friction. Markers are advected within the velocity field interpolated from the nodes. Marker properties are then passed back to the nodes to prepare the next solve of the conservation equations at the next time step. Markers also carry material stresses in order to solve Equation 32. In addition to being advected, these stress components are also rotated according to the local rotation rate determined from the velocity field at each timestep (Gerya, 2010).
# 4. Application to a 2-D Plane-Strain Numerical Press

## 4.1. Conservation Equations and Numerical Methodology

In order to perform 2-D simulations of material deformation governed by the SCAM constitutive equations, we adapt the long-term tectonic modeling code SiStER (Simple Stokes solver with Exotic Rheologies, Olive et al. (2016)), which solves for conservation of mass, momentum (and energy if needed), in a 2-D continuum assuming elastic incompressibility and planar deformation. Conservation of mass and momentum write:

$$ \frac{\partial v_i}{\partial x_i} = 0, \qquad (29) $$

and

$$ \frac{\partial s_{ij}}{\partial x_j} - \frac{\partial P}{\partial x_i} + \rho g_i = 0, \qquad (30) $$

where $v_i$ are velocities, $P = -1/3 \sigma_{kk}$ is pressure, $\rho$ is density and $g_i$ the gravitational acceleration. Deviatoric stresses $s_{ij}$ are related to velocities in Equation 30 using a Maxwell visco-elastic constitutive relationship between deviatoric stresses and strain rates (e.g., Gerya, 2010; Moresi et al., 2003):

$$ \dot{e}_{ij} = \frac{1}{2G} \dot{s}_{ij} + \frac{1}{2\eta} s_{ij} \qquad (31) $$

where $\dot{s}_{ij}$ is discretized using a first-order backward finite difference scheme with time step $\Delta t$, so that the deviatoric stress at time $t$ becomes

$$ s_{ij}^t = 2Z\eta\dot{e}_{ij}^t + (1-Z)s_{ij}^{t-\Delta t}, \qquad (32) $$

with $\dot{e}_{ij} = (\partial v_i / \partial x_j + \partial v_j / \partial x_i) / 2$ and

$$ Z = \frac{G\Delta t}{\eta + G\Delta t}. \qquad (33) $$

The effective viscosity in Equation 31 can represent a range of rheologies. A very high value sets a very long Maxwell time, which effectively renders the material elastic. In the viscous regime, $\eta$ can represent brittle plasticity (e.g., as detailed in Section 4.2.2), or a specific creep mechanism of known flow law. In practice, $\eta$ is constructed as the harmonic average of several viscosities, each representing individual flow mechanisms.

The mass Equation 29 and momentum Equation 30 conservation equations, expressed in terms of velocities Equation 32 and pressure, are discretized with a conservative finite difference scheme formulated on a staggered grid (e.g., Gerya & Yuen, 2003). This leads to a linear system that is solved for velocities and pressure over the entire domain using a direct solver. Retroactions between the viscosity and velocity fields require the use of non-linear iterations (here approximate-Newton, described as Algorithm 2 in Spiegelman et al. (2016)) to reach convergence, which is assessed by comparing the L2 norm of the residual vector to a specified tolerance (relative tolerance between $10^{-7}$ and $10^{-2}$, see readme documentation in the code repository linked in the Acknowledgments section).

Once a reasonably converged solution is found, the time evolution is performed explicitly using the time step $\Delta t$ introduced in Equation 32. This is specifically done by advecting Lagrangian markers which carry material properties such as density, viscosity and friction. Markers are advected within the velocity field interpolated from the nodes. Marker properties are then passed back to the nodes to prepare the next solve of the conservation equations at the next time step. Markers also carry material stresses in order to solve Equation 32. In addition to being advected, these stress components are also rotated according to the local rotation rate determined from the velocity field at each timestep (Gerya, 2010).
## 4.2. Numerical Implementation of the SCAM Rheology

The implementation of the SCAM model in the 2-D code was performed as described in the following sub-sections. First, a damage property and its evolution rules are implemented, along with the shear modulus dependence on damage. Then, a smooth transition to long term plastic behavior in fully damaged parts of the material is introduced.

### 4.2.1. Damage Growth and Viscosity

The damage state *D* is added as an additional variable discretized on both markers and nodes. Its evolution Equation 8 is solved with a finite difference method. The damage rate $\dot{D}$ and its associated viscosity $\eta_D$ (Equation 5) are evaluated on nodes at each non-linear iteration using previous stresses and interpolated damage values from markers. The shear modulus is also altered according to the damage state.

When stepping through time, marker damage is incremented by interpolating the damage rate from nodes to markers. Due to the non-linearity of Equation 8, damage is prone to catastrophic growth, which can be challenging for a numerical solver. We therefore adapt the time step to the dynamics of damage growth by limiting the maximum increment of damage on a node at each time iteration by an amount $\Delta D_{max}$.

### 4.2.2. Switching From Damaged to Plastic Rheology After Crack Coalescence

*D* values approaching 1 can be thought of as a state when the rock loses its macroscopic cohesion through crack coalescence. Our damage rheology is not well suited to represent large strains that may develop beyond this point, for example, within localized fault zones. Mohr-Coulomb plasticity, on the other hand, is perfectly relevant to model the frictional rheology of such fault gouges. Crack coalescence is however a necessary condition to the formation of macroscopic fault zones, such that damage growth up to 1 has to precede Mohr-Coulomb plastic deformation. Our damage model being formulated as an effective Maxwell rheology, we choose to retain this framework in our implementation of plasticity. We therefore implement a continuous effective viscosity that smoothly switches from our damage viscosity $\eta_D$ to the standard plastic viscosity $\eta_p$ (Duretz et al., 2021) as *D* approaches 1.

Because in our micromechanical model crack normals lie in the {σ₁, σ₃} plane, plastic deformation beyond coalescence should be confined to that same plane, and can therefore be modeled with a Mohr-Coulomb yield criterion ensuring that

$$s_{II} \le \sigma_y, \quad (34)$$

where $s_{II} = \sqrt{J_2(\sigma)}$, with $J_2(\sigma) = 1/2 s_{ij}s_{ij}$ the second invariant of the deviatoric stress tensor. $s_{II}$ is also the radius of Mohr's circle in 2-D incompressible plane strain:

$$s_{II} = \frac{1}{2}(\sigma_3 - \sigma_1). \quad (35)$$

In Equation 34, the plastic yield stress $\sigma_y$ writes:

$$\sigma_y = \sin\phi_m P + \cos\phi_m C_m. \quad (36)$$

The yield stress is a function of the macroscopic friction angle $\phi_m = \arctan \mu_m$ and cohesion $C_m$, as well as of the in-plane pressure $P$:

$$P = -\frac{1}{2}(\sigma_1 + \sigma_3), \quad (37)$$

which in elastically incompressible 2-D plane strain is also equal to total pressure. Satisfying the Mohr-Coulomb yield criterion within a Maxwell visco-elastic framework can be done through an effective “plastic viscosity” approach (e.g., Gerya (2010); Duretz et al. (2021)). As fully damaged areas become incompressible elastic-plastic zones, Equation 4 becomes:
## 4.2. Numerical Implementation of the SCAM Rheology

The implementation of the SCAM model in the 2-D code was performed as described in the following sub-sections. First, a damage property and its evolution rules are implemented, along with the shear modulus dependence on damage. Then, a smooth transition to long term plastic behavior in fully damaged parts of the material is introduced.

### 4.2.1. Damage Growth and Viscosity

The damage state $D$ is added as an additional variable discretized on both markers and nodes. Its evolution Equation 8 is solved with a finite difference method. The damage rate $\dot{D}$ and its associated viscosity $\eta_D$ (Equation 5) are evaluated on nodes at each non-linear iteration using previous stresses and interpolated damage values from markers. The shear modulus is also altered according to the damage state.

When stepping through time, marker damage is incremented by interpolating the damage rate from nodes to markers. Due to the non-linearity of Equation 8, damage is prone to catastrophic growth, which can be challenging for a numerical solver. We therefore adapt the time step to the dynamics of damage growth by limiting the maximum increment of damage on a node at each time iteration by an amount $\Delta D_{max}$.

### 4.2.2. Switching From Damaged to Plastic Rheology After Crack Coalescence

$D$ values approaching 1 can be thought of as a state when the rock loses its macroscopic cohesion through crack coalescence. Our damage rheology is not well suited to represent large strains that may develop beyond this point, for example, within localized fault zones. Mohr-Coulomb plasticity, on the other hand, is perfectly relevant to model the frictional rheology of such fault gouges. Crack coalescence is however a necessary condition to the formation of macroscopic fault zones, such that damage growth up to 1 has to precede Mohr-Coulomb plastic deformation. Our damage model being formulated as an effective Maxwell rheology, we choose to retain this framework in our implementation of plasticity. We therefore implement a continuous effective viscosity that smoothly switches from our damage viscosity $\eta_D$ to the standard plastic viscosity $\eta_p$ (Duretz et al., 2021) as $D$ approaches 1.

Because in our micromechanical model crack normals lie in the $\{\sigma_1, \sigma_3\}$ plane, plastic deformation beyond coalescence should be confined to that same plane, and can therefore be modeled with a Mohr-Coulomb yield criterion ensuring that

$$s_{II} \le \sigma_y, \qquad (34)$$

where $s_{II} = \sqrt{J_2(\sigma)}$, with $J_2(\sigma) = 1/2 s_{ij}s_{ij}$ the second invariant of the deviatoric stress tensor. $s_{II}$ is also the radius of Mohr's circle in 2-D incompressible plane strain:

$$s_{II} = \frac{1}{2}(\sigma_3 - \sigma_1). \qquad (35)$$

In Equation 34, the plastic yield stress $\sigma_y$ writes:

$$\sigma_y = \sin\phi_m P + \cos\phi_m C_m. \qquad (36)$$

The yield stress is a function of the macroscopic friction angle $\phi_m = \arctan \mu_m$ and cohesion $C_m$, as well as of the in-plane pressure $P$:

$$P = -\frac{1}{2}(\sigma_1 + \sigma_3), \qquad (37)$$

which in elastically incompressible 2-D plane strain is also equal to total pressure. Satisfying the Mohr-Coulomb yield criterion within a Maxwell visco-elastic framework can be done through an effective “plastic viscosity” approach (e.g., Gerya (2010); Duretz et al. (2021)). As fully damaged areas become incompressible elastic-plastic zones, Equation 4 becomes:
$$ \dot{e}_{ij} = \frac{\dot{s}_{ij}}{2G_0 f(D)} + \frac{s_{ij}}{2\eta_p}. \qquad (38) $$

The plastic viscosity $\eta_p$ is set to guarantee that stresses satisfy the Mohr-Coulomb criterion Equation 34 once “broken” material starts behaving as a plastic fault zone. If $s_{II}$ lies below $\sigma_y$, $\eta_p$ is effectively infinite (plasticity is not activated), otherwise $\eta_p$ reads (see Appendix B):

$$ \eta_p = \frac{s_{II}}{2\left(\dot{e}_{II} - \frac{\sigma_y - s_{II}}{2G_0 f(D) \Delta t}\right)}, \qquad (39) $$

where $\dot{e}_{II} = \sqrt{J_2(\dot{e})}$ is the second invariant of the strain rate tensor. A smooth transition from damage to plastic viscosity is implemented using a hyperbolic tangent function $S(x)$ that goes from 0 to 1 as its argument goes from negative to positive. We set the lowest viscosity to the plastic viscosity and thus write the continuously differentiable effective viscosity:

$$ \eta_{eff} = S(\eta_D - \eta_p)\eta_D + (1 - S(\eta_D - \eta_p))\eta_p. \qquad (40) $$

To ensure that the effective viscosity remains plastic when material is fully damaged, we set $\eta_D(D = 1)$ to the smallest viscosity that can be resolved by our numerical solver (see Section 4.3).

Including the transition to large-strain plasticity, the complete set of differential equations that constitute the SCAM model can be summarized as

$$ \dot{e}_{ij} = \frac{\dot{s}_{ij}}{2G_0 f(D)} + \frac{s_{ij}}{2\eta_{eff}}. \qquad (41) $$

When subjected to loading, the material first responds elastically with $D = D_i$ (its initial damage state $\ge D_0$), until $K_i$ becomes positive and damage starts growing (Figure 3). At that moment damage-driven alteration of the shear modulus generates an effective damage viscosity which affects the material behavior. Up to peak stress, the damage viscosity is greater than the plastic viscosity since it allows stress build-up, therefore $\eta_{eff} \sim \eta_D$. During the post-peak stress drop, $\eta_D$ quickly drops below the plastic viscosity which then becomes the effective viscosity. This allows the accumulation of large strains under stresses capped by the Mohr-Coulomb yield stress. Said yield stress is computed according to Equation 36. In the following, we adopt a macroscopic friction angle that matches the frictional properties of the shear defects, that is, $\phi_m = \tan^{-1}\mu$. We also assume that the fully damaged material is cohesionless, that is, $C_m = 0$. It should be noted that a fully damaged material may also return to an elastic behavior if $s_{II}$ happens to drop below $\sigma_y$. However, its shear modulus will have been permanently reduced by damage ($G = \gamma G_0$).

### 4.3. 2-D Setup: The Numerical Press

We construct a 2-D plane-strain analog to the triaxial experimental setup described in Section 3, following the geometry shown in Figure 9. This allows us to simulate constant strain rate deformation of a Westerly granite sample, with micromechanical properties determined in Section 3.3 (Table 1), up to large strains and including localization. The axial symmetry of triaxial tests allows us to only consider half of the sample's cross-section. Our geometry thus consists of a half-sample 10 cm tall and 2 cm wide on the right side of a wider box (10.5 × 6 cm) that includes confining fluid left of the sample, and a 0.5 mm-thick “piston” above the sample (Figure 9). Constant strain rate conditions are enforced by pushing material inward from the top of the domain at a constant velocity. The piston is here to ensure that new material flowing in during deformation is not of sample type. Outward velocities are prescribed along the left boundary to preserve a constant volume in the computational domain. The confining fluid is modeled as a low-viscosity Newtonian medium, with pressure imposed at the lower left corner of the domain. Gravity is ignored. The initial sample damage fluctuates spatially between $D_i = 0.136$ and 0.236 with an isotropic Perlin noise structure that represents material heterogeneities. The spatial domain is discretized using cell sizes of 2 × 0.5 mm within 3 cm of the left wall, and 0.5 × 0.5 mm within 3 cm of the right wall, that is, the part of the domain containing the sample.
$$ \dot{e}_{ij} = \frac{\dot{s}_{ij}}{2G_0 f(D)} + \frac{s_{ij}}{2\eta_p}. \qquad (38) $$

The plastic viscosity $\eta_p$ is set to guarantee that stresses satisfy the Mohr-Coulomb criterion Equation 34 once “broken” material starts behaving as a plastic fault zone. If $s_{II}$ lies below $\sigma_y$, $\eta_p$ is effectively infinite (plasticity is not activated), otherwise $\eta_p$ reads (see Appendix B):

$$ \eta_p = \frac{s_{II}}{2\left(\dot{e}_{II} - \frac{\sigma_y - s_{II}}{2G_0 f(D) \Delta t}\right)}, \qquad (39) $$

where $\dot{e}_{II} = \sqrt{J_2(\dot{e})}$ is the second invariant of the strain rate tensor. A smooth transition from damage to plastic viscosity is implemented using a hyperbolic tangent function $S(x)$ that goes from 0 to 1 as its argument goes from negative to positive. We set the lowest viscosity to the plastic viscosity and thus write the continuously differentiable effective viscosity:

$$ \eta_{eff} = S(\eta_D - \eta_p)\eta_D + (1 - S(\eta_D - \eta_p))\eta_p. \qquad (40) $$

To ensure that the effective viscosity remains plastic when material is fully damaged, we set $\eta_D(D = 1)$ to the smallest viscosity that can be resolved by our numerical solver (see Section 4.3).

Including the transition to large-strain plasticity, the complete set of differential equations that constitute the SCAM model can be summarized as

$$ \dot{e}_{ij} = \frac{\dot{s}_{ij}}{2G_0 f(D)} + \frac{s_{ij}}{2\eta_{eff}}. \qquad (41) $$

When subjected to loading, the material first responds elastically with $D = D_i$ (its initial damage state $\ge D_0$), until $K_i$ becomes positive and damage starts growing (Figure 3). At that moment damage-driven alteration of the shear modulus generates an effective damage viscosity which affects the material behavior. Up to peak stress, the damage viscosity is greater than the plastic viscosity since it allows stress build-up, therefore $\eta_{eff} \sim \eta_D$. During the post-peak stress drop, $\eta_D$ quickly drops below the plastic viscosity which then becomes the effective viscosity. This allows the accumulation of large strains under stresses capped by the Mohr-Coulomb yield stress. Said yield stress is computed according to Equation 36. In the following, we adopt a macroscopic friction angle that matches the frictional properties of the shear defects, that is, $\phi_m = \tan^{-1}\mu$. We also assume that the fully damaged material is cohesionless, that is, $C_m = 0$. It should be noted that a fully damaged material may also return to an elastic behavior if $s_{II}$ happens to drop below $\sigma_y$. However, its shear modulus will have been permanently reduced by damage ($G = \gamma G_0$).

### 4.3. 2-D Setup: The Numerical Press

We construct a 2-D plane-strain analog to the triaxial experimental setup described in Section 3, following the geometry shown in Figure 9. This allows us to simulate constant strain rate deformation of a Westerly granite sample, with micromechanical properties determined in Section 3.3 (Table 1), up to large strains and including localization. The axial symmetry of triaxial tests allows us to only consider half of the sample's cross-section. Our geometry thus consists of a half-sample 10 cm tall and 2 cm wide on the right side of a wider box (10.5 × 6 cm) that includes confining fluid left of the sample, and a 0.5 mm-thick “piston” above the sample (Figure 9). Constant strain rate conditions are enforced by pushing material inward from the top of the domain at a constant velocity. The piston is here to ensure that new material flowing in during deformation is not of sample type. Outward velocities are prescribed along the left boundary to preserve a constant volume in the computational domain. The confining fluid is modeled as a low-viscosity Newtonian medium, with pressure imposed at the lower left corner of the domain. Gravity is ignored. The initial sample damage fluctuates spatially between $D_i = 0.136$ and 0.236 with an isotropic Perlin noise structure that represents material heterogeneities. The spatial domain is discretized using cell sizes of 2 × 0.5 mm within 3 cm of the left wall, and 0.5 × 0.5 mm within 3 cm of the right wall, that is, the part of the domain containing the sample.
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Geometric setup used to simulate experimental deformation in 2D plane-strain. Exploiting the axial symmetry of triaxial experiments, we only model the left-half of a cross section containing the sample axis. The half-sample has a length of 10 cm and a width of 2 cm. It is initially seeded with a noisy damage field. A small piston one cell tall (with the same mechanical properties as the sample) is pushed in above the sample. Left of this assemblage lies the confining fluid kept at a constant pressure. Top and left borders are associated with Dirichlet boundary conditions on the velocity component normal to the borders, while the tangential components are left free. The right and bottom borders are free slip. Axial stress is evaluated by averaging vertical stresses in a 0.5 cm thick slice at the bottom of the sample.</figcaption></figure>

The convergence of Stokes solvers being very sensitive to viscosity contrasts, we restrict the maximum variation in viscosity across the domain to five orders of magnitude. We do so by setting an upper bound on viscosity $\eta_{max}$ such that its associated Maxwell time $\eta_{max}/G$ is 50 times longer than the time required to elastically reach the peak stress at the imposed axial strain rate. This large viscosity is initially assigned to the sample, rendering it effectively elastic at the beginning of the simulation. A lower bound on viscosity in the numerical domain is obtained through $\eta_{min} = 10^{-5} \times \eta_{max}$. This low viscosity is assigned to the confining fluid, ensuring that it behaves viscously throughout the simulation. Because the damage viscosity $\eta_D$ drops significantly as damage accumulates, the effective viscosity of the sample will decrease as it begins to fail. As $\eta_D$ approaches $\eta_p$, a smooth transition toward plastic viscosity is performed over a viscosity range $|\eta_D - \eta_p| \approx \eta_{min}/50$. Regardless of the viscosity transition, damage keeps increasing until reaching 1. At this point, it stops evolving and $\eta_D$ is fixed at $\eta_{min}$. Once parts of the sample are fully damaged, they effectively behave as a Mohr-Coulomb plastic solid with no cohesion and the same (macroscopic) friction coefficient as that determined to act on the microscopic shear defects (0.7).

Finally, the comparisons of simulated constant strain rate experiments with laboratory data requires the evaluation of macroscopic axial strains and deviatoric stresses. The axial strains are measured by tracking the displacement of the top boundary of the sample through time, and normalizing it by the initial size of the sample. Axial deviatoric stresses are obtained by averaging the vertical deviatoric stress $s_{ax}$ in a horizontal “stress gauge,” that is, a 0.5 cm thick band at the bottom of the numerical sample (Figure 9), excluding a cell size length near the left boundary, to avoid any influence from interpolations at the interface between sample and fluid.

## 4.4. Results

Figure 10 shows results from the numerical press performed under a constant axial strain rate of $10^{-5} \text{ s}^{-1}$ (Panels a and b) and $10^{-15} \text{ s}^{-1}$ (Panel b), and confining pressures of 30, 80, and 150 MPa. Figure 10c illustrates the patterns of damage growth and plastic strain for the simulation performed under 150 MPa of confining pressure (black lines in Panels a and b), on the left and right halves of each snapshot, respectively. The timing of each snapshot is indicated by the numbers on the stress-strain curve in Panel b. Up to snapshot 3, damage grows in a distributed fashion, which smoothes the initial heterogeneities. Damage increases homogeneously up to ~0.45 during that stage. Just prior to the peak stress, damage growth starts to localize close to the sample border, forming fast-growing damage bands at angles of ~30° with respect to the compression direction. They develop within a strain range of less than 0.01% (snapshots 5, 6, 7). Plastic strain is estimated by integrating the second invariant of the inelastic deviatoric strain rate tensor through time (once Mohr-Coulomb plasticity has been activated), and accumulates within fully damaged bands. Plastic shear banding first lags behind damage banding. Once a sample-scale damage band has grown, it effectively becomes a plastic shear band. This process begins as the axial stress drops abruptly (snapshots 6–7). Interestingly, off-band distributed damage does not evolve significantly during shear band development.

Figure 10a compares the stress-strain curves produced in our 2-D simulations to the experiments of Wawersik and Brace (1971) described in Section 3.3, performed under the same conditions. To facilitate comparisons between a 2-D plane strain and an axisymmetric setup, we normalize the differential stress $\Delta\sigma = \sigma_3 - \sigma_1$ by its value at $K_I = 0$ and $D = D_0$, which is the criterion for the onset of tensile crack growth (even though the crack growth rate is infinitely slow at $K_I = O^+$):
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Geometric setup used to simulate experimental deformation in 2D plane-strain. Exploiting the axial symmetry of triaxial experiments, we only model the left-half of a cross section containing the sample axis. The half-sample has a length of 10 cm and a width of 2 cm. It is initially seeded with a noisy damage field. A small piston one cell tall (with the same mechanical properties as the sample) is pushed in above the sample. Left of this assemblage lies the confining fluid kept at a constant pressure. Top and left borders are associated with Dirichlet boundary conditions on the velocity component normal to the borders, while the tangential components are left free. The right and bottom borders are free slip. Axial stress is evaluated by averaging vertical stresses in a 0.5 cm thick slice at the bottom of the sample.</figcaption></figure>

The convergence of Stokes solvers being very sensitive to viscosity contrasts, we restrict the maximum variation in viscosity across the domain to five orders of magnitude. We do so by setting an upper bound on viscosity $\eta_{max}$ such that its associated Maxwell time $\eta_{max}/G$ is 50 times longer than the time required to elastically reach the peak stress at the imposed axial strain rate. This large viscosity is initially assigned to the sample, rendering it effectively elastic at the beginning of the simulation. A lower bound on viscosity in the numerical domain is obtained through $\eta_{min} = 10^{-5} \times \eta_{max}$. This low viscosity is assigned to the confining fluid, ensuring that it behaves viscously throughout the simulation. Because the damage viscosity $\eta_D$ drops significantly as damage accumulates, the effective viscosity of the sample will decrease as it begins to fail. As $\eta_D$ approaches $\eta_p$, a smooth transition toward plastic viscosity is performed over a viscosity range $|\eta_D - \eta_p| \approx \eta_{min}/50$. Regardless of the viscosity transition, damage keeps increasing until reaching 1. At this point, it stops evolving and $\eta_D$ is fixed at $\eta_{min}$. Once parts of the sample are fully damaged, they effectively behave as a Mohr-Coulomb plastic solid with no cohesion and the same (macroscopic) friction coefficient as that determined to act on the microscopic shear defects (0.7).

Finally, the comparisons of simulated constant strain rate experiments with laboratory data requires the evaluation of macroscopic axial strains and deviatoric stresses. The axial strains are measured by tracking the displacement of the top boundary of the sample through time, and normalizing it by the initial size of the sample. Axial deviatoric stresses are obtained by averaging the vertical deviatoric stress $s_{ax}$ in a horizontal “stress gauge,” that is, a 0.5 cm thick band at the bottom of the numerical sample (Figure 9), excluding a cell size length near the left boundary, to avoid any influence from interpolations at the interface between sample and fluid.

## 4.4. Results

Figure 10 shows results from the numerical press performed under a constant axial strain rate of $10^{-5} \text{ s}^{-1}$ (Panels a and b) and $10^{-15} \text{ s}^{-1}$ (Panel b), and confining pressures of 30, 80, and 150 MPa. Figure 10c illustrates the patterns of damage growth and plastic strain for the simulation performed under 150 MPa of confining pressure (black lines in Panels a and b), on the left and right halves of each snapshot, respectively. The timing of each snapshot is indicated by the numbers on the stress-strain curve in Panel b. Up to snapshot 3, damage grows in a distributed fashion, which smoothes the initial heterogeneities. Damage increases homogeneously up to ~0.45 during that stage. Just prior to the peak stress, damage growth starts to localize close to the sample border, forming fast-growing damage bands at angles of ~30° with respect to the compression direction. They develop within a strain range of less than 0.01% (snapshots 5, 6, 7). Plastic strain is estimated by integrating the second invariant of the inelastic deviatoric strain rate tensor through time (once Mohr-Coulomb plasticity has been activated), and accumulates within fully damaged bands. Plastic shear banding first lags behind damage banding. Once a sample-scale damage band has grown, it effectively becomes a plastic shear band. This process begins as the axial stress drops abruptly (snapshots 6–7). Interestingly, off-band distributed damage does not evolve significantly during shear band development.

Figure 10a compares the stress-strain curves produced in our 2-D simulations to the experiments of Wawersik and Brace (1971) described in Section 3.3, performed under the same conditions. To facilitate comparisons between a 2-D plane strain and an axisymmetric setup, we normalize the differential stress $\Delta\sigma = \sigma_3 - \sigma_1$ by its value at $K_I = 0$ and $D = D_0$, which is the criterion for the onset of tensile crack growth (even though the crack growth rate is infinitely slow at $K_I = O^+$):
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Results of the SCAM model simulations in a 2-D plane-strain setup using Westerly granite parameters under constant strain rate. The numerical samples are initialized with a noisy $D_i$ field. (a) Simulated stress-strain curves at $\dot{e}_{ax} = -10^{-5} s^{-1}$ and confining pressures of 30 (red), 80 (blue) and 150 MPa (black curve), to be compared with the corresponding experimental data (dots). (b) Simulation performed at $p_c = 150$ MPa described above compared to a simulation at a “tectonic” strain rate of $10^{-15} s^{-1}$. Differential stresses and deviatoric axial strains are normalized by their value when $K_f = 0$ at $D = D_0$ to allow comparison between triaxial experimental data and 2-D plane strain deformation (see Text). The red numbers in Panel b correspond to the damage (left) and plastic strain (right) snapshots shown in Panel (c).</figcaption></figure>

$$ \Delta\sigma^c = p_c \left( 1 - \frac{\sqrt{1 + \mu^2} + \mu}{\sqrt{1 + \mu^2} - \mu} \right) \qquad (42) $$

The above expression is obtained by applying the Mohr-Coulomb criterion to optimally oriented planes in a principal stress field. The axial deviatoric strain is then normalized by the deviatoric strain needed to reach $\Delta\sigma^c$ elastically with the reference shear modulus $G_0$ that corresponds to $D = D_0$:

$$ e_{ax}^c = \frac{1}{2G_0} s_{ax}^c \qquad (43) $$

In Equation 43, $s_{ax}^c$ is the deviatoric axial stress at the onset of microcraking, and is equal to $\Delta\sigma^c/2$ in a 2-D plane strain configuration, and to $(2/3)\Delta\sigma^c$ for a triaxial configuration. This non-dimensionalization of stresses and strains accounts for the fact that the mean stress, which impacts damage growth and the position of the peak stress, has a different expression in a triaxial versus plane-strain geometry. It ensures that experiments conducted with
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Results of the SCAM model simulations in a 2-D plane-strain setup using Westerly granite parameters under constant strain rate. The numerical samples are initialized with a noisy $D_i$ field. (a) Simulated stress-strain curves at $\dot{e}_{ax} = -10^{-5} s^{-1}$ and confining pressures of 30 (red), 80 (blue) and 150 MPa (black curve), to be compared with the corresponding experimental data (dots). (b) Simulation performed at $p_c = 150$ MPa described above compared to a simulation at a “tectonic” strain rate of $10^{-15} s^{-1}$. Differential stresses and deviatoric axial strains are normalized by their value when $K_f = 0$ at $D = D_0$ to allow comparison between triaxial experimental data and 2-D plane strain deformation (see Text). The red numbers in Panel b correspond to the damage (left) and plastic strain (right) snapshots shown in Panel (c).</figcaption></figure>

$$ \Delta\sigma^c = p_c \left( 1 - \frac{\sqrt{1 + \mu^2} + \mu}{\sqrt{1 + \mu^2} - \mu} \right) \qquad (42) $$

The above expression is obtained by applying the Mohr-Coulomb criterion to optimally oriented planes in a principal stress field. The axial deviatoric strain is then normalized by the deviatoric strain needed to reach $\Delta\sigma^c$ elastically with the reference shear modulus $G_0$ that corresponds to $D = D_0$:

$$ \epsilon_{ax}^c = \frac{1}{2G_0} s_{ax}^c \qquad (43) $$

In Equation 43, $s_{ax}^c$ is the deviatoric axial stress at the onset of microcraking, and is equal to $\Delta\sigma^c/2$ in a 2-D plane strain configuration, and to $(2/3)\Delta\sigma^c$ for a triaxial configuration. This non-dimensionalization of stresses and strains accounts for the fact that the mean stress, which impacts damage growth and the position of the peak stress, has a different expression in a triaxial versus plane-strain geometry. It ensures that experiments conducted with
the same parameters in either geometry will show the same non-dimensional tangent modulus and peak differential stress.

Our 2-D simulations are in good agreement with the experimental data from Wawersik and Brace (1971) (Figure 10a) at all three confining pressures. This shows that the parameters determined by fitting 0-D (point-wise) simulations to triaxial data produce sensible behavior when implemented in a 2-D “spatialized” geometry. Figure 10b shows our reference simulation at $\dot{e}_{ax} = 10^{-5} \text{ s}^{-1}$ and $p_c = 150 \text{ MPa}$ in black, compared to a simulation performed at a “tectonic” strain rate of $10^{-15} \text{ s}^{-1}$. The peak stress of the slower simulation is significantly lower than that of the reference simulation, with a loss of strength during macroscopic failure that is approximately divided by 2.

# 5. Discussion: A Brittle Constitutive Law Rooted in Micromechanics

## 5.1. Features of Brittle Deformation Captured by the SCAM Model

As illustrated in Sections 3 and 4, the SCAM model captures a range of features typical of brittle deformation revealed by laboratory experiments (Figure 1). These include: (a) the co-existence of several measures of rock strength, such as the intact strength and the residual (i.e., “pre-cut”) frictional strength, all of which depend on confining pressure (e.g., Byerlee (1978)); (b) the permanent weakening of elastic properties occurring prior to the peak stress; (c) the strain rate-dependence of brittle strength, which enables (d) the occurrence of brittle creep under constant imposed stress. Here we discuss the parameters of the SCAM framework that control these various macroscopic properties.

### 5.1.1. Microscopic Versus Macroscopic Strength, Elastic Weakening, and Strain Rate Dependence

The SCAM framework involves several thresholds of inelastic deformation. The first is when slip on pre-existing, small-scale shear defects becomes able to wedge open tensile wing cracks. It corresponds to $K_I = 0$ (Figures 3 and 6), and is closely related to the Mohr-Coulomb criterion, in the sense that opening wing cracks requires a greater differential stress under greater confining pressure (Costin, 1985). The second threshold is when cracks have sufficiently lengthened to transition from a non-interacting to an interacting regime. This aspect will be further detailed in Section 5.1.2. The third threshold is when $D$ reaches its maximum value $\sim 1$, at which point cracks coalesce into a macroscopic fault, which is modeled as a shear band with a macroscopic “bulk” friction equal to that acting on the microscopic shear defects, and no cohesion. In practice, the second and third thresholds occur in very close succession because the damage growth rate accelerates catastrophically as soon as cracks enter the interacting regime, especially under constant axial strain rate. Some amount of deformation is still required for the shear band to reach its steady-state stress after the third threshold (e.g., after snapshot 7 in Figure 10). After that, the macroscopically broken rock has a “residual” strength that is entirely set by its friction coefficient.

The material's elastic properties are altered by damage, causing pre-peak softening of the rock and permanent weakening of the shear modulus. The ratio of the fully damaged ($V_s^D = V_s(D=1)$) over the reference ($V_s^0 = V_s(D=D_0)$) shear-wave velocity can be related to the shear modulus weakening ratio $\gamma$, assuming small density variations, as: $V_s^D/V_s^0 \approx \sqrt{\gamma}$. Using $\gamma$ values inverted from Westerly granite and Darley Dale sandstone (Table 1) we obtain shear-wave velocity reductions of 34% and 47%, consistent with values measured in the damage zone of natural faults, which range from 20% to 50% (e.g., Karabulut & Bouchon, 2007; C. Wu et al., 2009), as well as laboratory tomography on granite showing a reduction in P-wave velocity of around 50% (Aben et al., 2019).

The SCAM model accounts for the temporal dependence of brittle deformation via a sub-critical crack growth law, which allows cracks to grow below the fracture toughness of the material (Costin, 1983, 1985). This assumption introduces a characteristic crack growth time $\alpha/\dot{\alpha}_0$ that is modulated by the stress intensity factor ($K_I$) and the Charles law exponent $n$ (Equation 7). These parameters themselves depend on ambient conditions such as moisture levels (Atkinson, 1979; Eppes & Keanini, 2017) or temperature (Heap et al., 2009). To first order, the strain rate dependence of the SCAM flow law reflects the ratio of the characteristic duration of the deformation of interest to the characteristic crack growth time. A very slow (“tectonic”) experiment will for example, leave ample time for cracks to grow, weaken the material and cause macroscopic failure, preventing the build-up of very large stresses (Figure 10b). Conversely, experiments conducted under laboratory strain rates will reveal greater peak
the same parameters in either geometry will show the same non-dimensional tangent modulus and peak differential stress.

Our 2-D simulations are in good agreement with the experimental data from Wawersik and Brace (1971) (Figure 10a) at all three confining pressures. This shows that the parameters determined by fitting 0-D (point-wise) simulations to triaxial data produce sensible behavior when implemented in a 2-D “spatialized” geometry. Figure 10b shows our reference simulation at $\dot{e}_{ax} = 10^{-5}$ s$^{-1}$ and $p_c = 150$ MPa in black, compared to a simulation performed at a “tectonic” strain rate of $10^{-15}$ s$^{-1}$. The peak stress of the slower simulation is significantly lower than that of the reference simulation, with a loss of strength during macroscopic failure that is approximately divided by 2.

# 5. Discussion: A Brittle Constitutive Law Rooted in Micromechanics

## 5.1. Features of Brittle Deformation Captured by the SCAM Model

As illustrated in Sections 3 and 4, the SCAM model captures a range of features typical of brittle deformation revealed by laboratory experiments (Figure 1). These include: (a) the co-existence of several measures of rock strength, such as the intact strength and the residual (i.e., “pre-cut”) frictional strength, all of which depend on confining pressure (e.g., Byerlee (1978)); (b) the permanent weakening of elastic properties occurring prior to the peak stress; (c) the strain rate-dependence of brittle strength, which enables (d) the occurrence of brittle creep under constant imposed stress. Here we discuss the parameters of the SCAM framework that control these various macroscopic properties.

### 5.1.1. Microscopic Versus Macroscopic Strength, Elastic Weakening, and Strain Rate Dependence

The SCAM framework involves several thresholds of inelastic deformation. The first is when slip on pre-existing, small-scale shear defects becomes able to wedge open tensile wing cracks. It corresponds to $K_I = 0$ (Figures 3 and 6), and is closely related to the Mohr-Coulomb criterion, in the sense that opening wing cracks requires a greater differential stress under greater confining pressure (Costin, 1985). The second threshold is when cracks have sufficiently lengthened to transition from a non-interacting to an interacting regime. This aspect will be further detailed in Section 5.1.2. The third threshold is when $D$ reaches its maximum value $\sim 1$, at which point cracks coalesce into a macroscopic fault, which is modeled as a shear band with a macroscopic “bulk” friction equal to that acting on the microscopic shear defects, and no cohesion. In practice, the second and third thresholds occur in very close succession because the damage growth rate accelerates catastrophically as soon as cracks enter the interacting regime, especially under constant axial strain rate. Some amount of deformation is still required for the shear band to reach its steady-state stress after the third threshold (e.g., after snapshot 7 in Figure 10). After that, the macroscopically broken rock has a “residual” strength that is entirely set by its friction coefficient.

The material's elastic properties are altered by damage, causing pre-peak softening of the rock and permanent weakening of the shear modulus. The ratio of the fully damaged ($V_s^D = V_s(D=1)$) over the reference ($V_s^0 = V_s(D=D_0)$) shear-wave velocity can be related to the shear modulus weakening ratio $\gamma$, assuming small density variations, as: $V_s^D/V_s^0 \approx \sqrt{\gamma}$. Using $\gamma$ values inverted from Westerly granite and Darley Dale sandstone (Table 1) we obtain shear-wave velocity reductions of 34% and 47%, consistent with values measured in the damage zone of natural faults, which range from 20% to 50% (e.g., Karabulut & Bouchon, 2007; C. Wu et al., 2009), as well as laboratory tomography on granite showing a reduction in P-wave velocity of around 50% (Aben et al., 2019).

The SCAM model accounts for the temporal dependence of brittle deformation via a sub-critical crack growth law, which allows cracks to grow below the fracture toughness of the material (Costin, 1983, 1985). This assumption introduces a characteristic crack growth time $\alpha/\dot{\alpha}_0$ that is modulated by the stress intensity factor ($K_I$) and the Charles law exponent $n$ (Equation 7). These parameters themselves depend on ambient conditions such as moisture levels (Atkinson, 1979; Eppes & Keanini, 2017) or temperature (Heap et al., 2009). To first order, the strain rate dependence of the SCAM flow law reflects the ratio of the characteristic duration of the deformation of interest to the characteristic crack growth time. A very slow (“tectonic”) experiment will for example, leave ample time for cracks to grow, weaken the material and cause macroscopic failure, preventing the build-up of very large stresses (Figure 10b). Conversely, experiments conducted under laboratory strain rates will reveal greater peak
stress (e.g., Figure 7 of Costin (1983)). Interestingly, the strain rate dependence of brittle deformation implies that wide portions of the upper crust should behave in an effectively viscous fashion (with viscosity η_D) when undergoing progressive failure (i.e., prior to localization). This behavior must however be inherently transient because the amount of damage a rock can withstand before macroscopic failure is necessarily finite. A competition between crack growth and crack healing processes may prolong this distributed viscous deformation phase, but is beyond the scope of the present study.

### 5.1.2. Retroactions Between Damage and Damage Growth Rate

The transition between a regime of lengthening but non-interacting wing cracks, and one of catastrophically interacting long cracks (Figure 3) is at the heart of many macroscopic behaviors manifested by the SCAM model. As detailed in Section 2.3, the stress intensity factor ($K_I$) at the tip of wing cracks is constructed as the sum of three terms (Equation 9). The first two terms lead to a decrease in $K_I$ as wing cracks lengthen, that is, as D increases. This corresponds to the isolated crack regime. The third term has the opposite effect: increasing D increases $K_I$, and thus the damage growth rate through Charles' law (Equation 7). This last term becomes dominant at larger values of D and describes the interacting crack regime. The transition between successive regimes is closely related to the convexity of $K_I$ as a function of D, and its dependence on the evolving differential stress, as illustrated in Figures 3e–3g.

In constant strain rate experiments, damage starts to grow when $K_I$ becomes positive. This is made possible by elastic loading raising the differential stress at constant initial damage state $D_I$ (vertical trajectories in Figure 6). When an experiment is started with a low $D_I$, damage growth first occurs in the non-interacting regime, where $\partial K_I / \partial D < 0$ (e.g., purple trajectory in Figure 6). Damage growth rates are initially very slow, because $K_I / K_{IC}$ raised to a large exponent (n in Equation 7) gives an extremely slow crack growth speed when $K_I$ barely exceeds 0. The damage viscosity is initially very high, and the material continues to behave elastically. As damage increases, both the shear modulus and the damage viscosity decrease because of the decreasing f(D) and $f^2(D)$ terms in Equations 2 and 5. This leads to pre-peak softening of the stress-strain curve. The crack growth rate—strongly controlled by $K_I$—is the sole mechanism that can lead to a stress rate decrease. It thus competes with the elastic stress rate increase imparted by far-field loading, controlling the stress level at which material softening occurs. In Figure 6, this manifests as trajectories aligning on a contour of constant $K_I$, which is greater for a greater imposed strain rate.

Then, the system transitions to the interacting regime where $\partial K_I / \partial D > 0$. The regime transition as illustrated in Figure 6a connects all the differential stress maxima spanning all values of $K_I$ between 0 and $K_{IC}$. We write $D_c$ the “critical” damage value that marks this regime transition. $D_c$ decreases with increasing differential stress (vertical black dashed curve in Figure 6a), and verifies:

$$ \left. \frac{\partial K_I(D)}{\partial D} \right|_{D=D_c(\Delta\sigma)} = 0. \qquad (44) $$

Because prior to crossing the regime transition stress trajectories align close to an iso-$K_I$ contour (which depends on strain rate), the value of $D_c$ can be thought of as a decreasing function of strain rate. $D_c$ is bounded by the value of damage that maximizes the differential stress at $K_I = K_{IC}$ (here ~0.42), and the value that maximizes stress at $K_I = 0$ (here ~0.5). These end-member cases respectively represent a very fast strain rate experiment, in which cracks would grow critically (at elastic wave speeds), and an extremely long and slow experiment in which cracks can grow sub-critically at $K_I \sim 0$. The upper and lower bound on $D_c$ happen to be close to each other, yielding a narrow range of critical damage values (~0.42 – 0.5 in Figure 6a). When damage exceeds $D_c$, the system enters the interacting regime, in which an increase in D increases $\dot{D}$ at constant stress, thereby accelerating the reduction of the shear modulus (Equation 2) and damage viscosity $\eta_D$ (Equation 5). The material can no longer accumulate stress, and stresses decrease below their peak value. At this point the stress trajectories in Figure 6a begin to deviate from an iso-$K_I$. For our best-fitting set of parameters, $K_I$ increases drastically as D exceeds $D_c$, which manifests as a sharp stress drop as D approaches 1 (Figure 6).

An interesting consequence of the fact that pre-peak stress trajectories tend to first align on the same iso-$K_I$ contour regardless of initial damage state is that they all experience a regime transition at the same $D_c$ and at the same peak differential stress (for a given imposed strain rate). In Figure 6b, this manifests as peak stress
stress (e.g., Figure 7 of Costin (1983)). Interestingly, the strain rate dependence of brittle deformation implies that wide portions of the upper crust should behave in an effectively viscous fashion (with viscosity η_D) when undergoing progressive failure (i.e., prior to localization). This behavior must however be inherently transient because the amount of damage a rock can withstand before macroscopic failure is necessarily finite. A competition between crack growth and crack healing processes may prolong this distributed viscous deformation phase, but is beyond the scope of the present study.

### 5.1.2. Retroactions Between Damage and Damage Growth Rate

The transition between a regime of lengthening but non-interacting wing cracks, and one of catastrophically interacting long cracks (Figure 3) is at the heart of many macroscopic behaviors manifested by the SCAM model. As detailed in Section 2.3, the stress intensity factor (K_I) at the tip of wing cracks is constructed as the sum of three terms (Equation 9). The first two terms lead to a decrease in K_I as wing cracks lengthen, that is, as D increases. This corresponds to the isolated crack regime. The third term has the opposite effect: increasing D increases K_I, and thus the damage growth rate through Charles' law (Equation 7). This last term becomes dominant at larger values of D and describes the interacting crack regime. The transition between successive regimes is closely related to the convexity of K_I as a function of D, and its dependence on the evolving differential stress, as illustrated in Figures 3e–3g.

In constant strain rate experiments, damage starts to grow when K_I becomes positive. This is made possible by elastic loading raising the differential stress at constant initial damage state D_I (vertical trajectories in Figure 6). When an experiment is started with a low D_I, damage growth first occurs in the non-interacting regime, where ∂K_I/∂D < 0 (e.g., purple trajectory in Figure 6). Damage growth rates are initially very slow, because K_I/K_{IC} raised to a large exponent (n in Equation 7) gives an extremely slow crack growth speed when K_I barely exceeds 0. The damage viscosity is initially very high, and the material continues to behave elastically. As damage increases, both the shear modulus and the damage viscosity decrease because of the decreasing f(D) and f^2(D) terms in Equations 2 and 5. This leads to pre-peak softening of the stress-strain curve. The crack growth rate—strongly controlled by K_I—is the sole mechanism that can lead to a stress rate decrease. It thus competes with the elastic stress rate increase imparted by far-field loading, controlling the stress level at which material softening occurs. In Figure 6, this manifests as trajectories aligning on a contour of constant K_I, which is greater for a greater imposed strain rate.

Then, the system transitions to the interacting regime where ∂K_I/∂D > 0. The regime transition as illustrated in Figure 6a connects all the differential stress maxima spanning all values of K_I between 0 and K_{IC}. We write D_c the “critical” damage value that marks this regime transition. D_c decreases with increasing differential stress (vertical black dashed curve in Figure 6a), and verifies:

$$ \left. \frac{\partial K_I(D)}{\partial D} \right|_{D=D_c(\Delta\sigma)} = 0. \qquad (44) $$

Because prior to crossing the regime transition stress trajectories align close to an iso-K_I contour (which depends on strain rate), the value of D_c can be thought of as a decreasing function of strain rate. D_c is bounded by the value of damage that maximizes the differential stress at K_I = K_{IC} (here ~0.42), and the value that maximizes stress at K_I = 0 (here ~0.5). These end-member cases respectively represent a very fast strain rate experiment, in which cracks would grow critically (at elastic wave speeds), and an extremely long and slow experiment in which cracks can grow sub-critically at K_I ~ 0. The upper and lower bound on D_c happen to be close to each other, yielding a narrow range of critical damage values (~0.42 – 0.5 in Figure 6a). When damage exceeds D_c, the system enters the interacting regime, in which an increase in D increases Ddot at constant stress, thereby accelerating the reduction of the shear modulus (Equation 2) and damage viscosity η_D (Equation 5). The material can no longer accumulate stress, and stresses decrease below their peak value. At this point the stress trajectories in Figure 6a begin to deviate from an iso-K_I. For our best-fitting set of parameters, K_I increases drastically as D exceeds D_c, which manifests as a sharp stress drop as D approaches 1 (Figure 6).

An interesting consequence of the fact that pre-peak stress trajectories tend to first align on the same iso-K_I contour regardless of initial damage state is that they all experience a regime transition at the same D_c and at the same peak differential stress (for a given imposed strain rate). In Figure 6b, this manifests as peak stress
magnitudes that are largely insensitive to any value of initial damage lower than $D_c$. This feature of the model is consistent with experimental results from Wang et al. (2013), who found similar peak strength in samples initially subjected to varying degrees of thermal cracking, which we interpret to represent varying $D_I$ (i.e., varying wing crack lengths at fixed shear defect size and density). On the other hand, if an experiment is started with a damage state that exceeds $D_c$, the system will entirely bypass the non-interacting regime and will display very little post-peak softening (e.g., light orange trajectories in Figure 6a). In this case, the degree of initial damage affects the position of the peak stress.

We note that in a few instances, large values of damage do not lead to a catastrophic stress drop. This occurs for example, for low values of $\gamma$ (Figure 5h), or a high value of $D_0$ (Figure 5f), where stresses slowly decay over a few percent of axial strain. The only way to prevent a catastrophic stress decrease is for $K_I$ to decrease as $D$ approaches 1. In Figure 6a, this would manifest as a steeply decreasing stress trajectory that crosses iso-$K_I$ contours for $D > D_c$. The slope of a stress trajectory in $(D, \Delta\sigma)$ space is equal to $\frac{\partial\Delta\sigma}{\partial D} = \frac{3}{2} \dot{s}_{ax}$. Post peak, the deviatoric axial stress rate (Equation 27) is increasingly dominated by the viscous term, as $\dot{D}$ accelerates. Neglecting the elastic term in Equation 27 yields:

$$ \frac{\partial \Delta \sigma}{\partial D} = -\frac{3}{2} \frac{|f'(D)|}{f(D)} s_{ax}. \qquad (45) $$

Using the equations for $f'(D)$ Equation 24 and $f(D)$ Equation 23, it can be seen that a low value of $\gamma$ leads to a very steep $\frac{\partial\Delta\sigma}{\partial D}$. This likely explains the gentler stress drops shown in Figure 5h. We suspect that a large $D_0$ leads to a similar effect on $\frac{\partial\Delta\sigma}{\partial D}$ and accounts for the progressive stress drop in Figure 5f. It is noteworthy that our inversion for Westerly granite predicts a sharp stress drop, even though it relies on data that does not span the interacting crack regime ($D > D_c$, i.e., only pre-peak data is used in Figure 7a, and minimum creep strain rates in Figure 7b).

Our 2-D simulations help us assess how the change of crack growth regime affects the spatial pattern of damage as it transitions from distributed to localized (Figure 10c). In particular, the initial phase of distributed damage growth appears to coincide with the non-interacting regime. Damage increases uniformly, smoothing any pre-existing initial damage heterogeneity, and reaches a near constant value (~0.45) in the bulk rock when damage bands begin forming. We interpret this uniform value as related to $D_c$. Specifically, the distributed build-up of damage (snapshots 1, 2, and 3 in Figure 10c) proceeds in the non-interacting regime, in which damage is uniformly capped at ~$D_c$. Stress concentrations due to numerical noise or prescribed heterogeneities can however trigger the switch to the interacting regime in some portions of the sample (along the sides in snapshot 5 of Figure 10), leading to the localization of damage bands. Interestingly, the stable, uniform growth of damage when $D < D_c$ is probably the reason for the good agreement between our 2-D simulations and our 0-D models, which are by definition “homogeneous” (Figure 10a). The post-peak behavior predicted by the SCAM model is however significantly different in 2-D versus 0-D simulations. This is because it is driven by retro-actions between damage localization within a band and the stress field of the surrounding rock that cannot be captured in a pointwise model.

The distinct regimes of crack growth are also responsible for the two stages of creep observed in our constant stress simulations (Figures 4c and 4d, see Section 3), as well as in brittle creep experiments (Figure 1). In the representation of Figure 6a, a constant stress experiment simply maps as a horizontal line starting from any point of the constant strain rate trajectories prior to the peak stress. In order to break the material under constant stress, the $K_I$ trajectory has to remain in the domain $K_I > 0$ up to $D = 1$. There thus exists a threshold in differential stress $\Delta\sigma_{bc}$ that must be met for the sample to fail macroscopically. Otherwise, the accumulation of damage under constant stress will decrease $K_I$ all the way to negative values, inhibiting the growth of further damage. This threshold corresponds to the largest differential stress able to produce a stress concentration factor equal to zero. It is represented in Figure 6a by the summit of the dashed contour of $K_I = 0$, and corresponds to a value of around 720 MPa for a confining pressure of 150 MPa with our inverted Westerly granite parameters. Overall, $\Delta\sigma_{bc}$ can be thought of as a theoretical minimum strength of the rock, which is a function of confining pressure only (continuous blue line in Figure 11). If a brittle creep test is carried out under a constant differential stress above $\Delta\sigma_{bc}$, the damage state will eventually reach $D_c$ in a finite amount of time, then transition to the interacting regime that allows failure. In this case, the creep test will begin by a decrease in $K_I$ that manifests as a decrease in the macroscopic strain rate referred to as decelerating or primary creep (Figure 4c). We note that the minimum brittle
magnitudes that are largely insensitive to any value of initial damage lower than $D_c$. This feature of the model is consistent with experimental results from Wang et al. (2013), who found similar peak strength in samples initially subjected to varying degrees of thermal cracking, which we interpret to represent varying $D_I$ (i.e., varying wing crack lengths at fixed shear defect size and density). On the other hand, if an experiment is started with a damage state that exceeds $D_c$, the system will entirely bypass the non-interacting regime and will display very little post-peak softening (e.g., light orange trajectories in Figure 6a). In this case, the degree of initial damage affects the position of the peak stress.

We note that in a few instances, large values of damage do not lead to a catastrophic stress drop. This occurs for example, for low values of $\gamma$ (Figure 5h), or a high value of $D_0$ (Figure 5f), where stresses slowly decay over a few percent of axial strain. The only way to prevent a catastrophic stress decrease is for $K_I$ to decrease as $D$ approaches 1. In Figure 6a, this would manifest as a steeply decreasing stress trajectory that crosses iso-$K_I$ contours for $D > D_c$. The slope of a stress trajectory in $(D, \Delta\sigma)$ space is equal to $\frac{\partial\Delta\sigma}{\partial D} = \frac{3}{2} \dot{s}_{ax}$. Post peak, the deviatoric axial stress rate (Equation 27) is increasingly dominated by the viscous term, as $\dot{D}$ accelerates. Neglecting the elastic term in Equation 27 yields:

$$ \frac{\partial \Delta \sigma}{\partial D} = -\frac{3}{2} \frac{|f'(D)|}{f(D)} s_{ax}. \quad (45) $$

Using the equations for $f'(D)$ Equation 24 and $f(D)$ Equation 23, it can be seen that a low value of $\gamma$ leads to a very steep $\frac{\partial\Delta\sigma}{\partial D}$. This likely explains the gentler stress drops shown in Figure 5h. We suspect that a large $D_0$ leads to a similar effect on $\frac{\partial\Delta\sigma}{\partial D}$ and accounts for the progressive stress drop in Figure 5f. It is noteworthy that our inversion for Westerly granite predicts a sharp stress drop, even though it relies on data that does not span the interacting crack regime ($D > D_c$, i.e., only pre-peak data is used in Figure 7a, and minimum creep strain rates in Figure 7b).

Our 2-D simulations help us assess how the change of crack growth regime affects the spatial pattern of damage as it transitions from distributed to localized (Figure 10c). In particular, the initial phase of distributed damage growth appears to coincide with the non-interacting regime. Damage increases uniformly, smoothing any pre-existing initial damage heterogeneity, and reaches a near constant value (~0.45) in the bulk rock when damage bands begin forming. We interpret this uniform value as related to $D_c$. Specifically, the distributed build-up of damage (snapshots 1, 2, and 3 in Figure 10c) proceeds in the non-interacting regime, in which damage is uniformly capped at ~$D_c$. Stress concentrations due to numerical noise or prescribed heterogeneities can however trigger the switch to the interacting regime in some portions of the sample (along the sides in snapshot 5 of Figure 10), leading to the localization of damage bands. Interestingly, the stable, uniform growth of damage when $D < D_c$ is probably the reason for the good agreement between our 2-D simulations and our 0-D models, which are by definition “homogeneous” (Figure 10a). The post-peak behavior predicted by the SCAM model is however significantly different in 2-D versus 0-D simulations. This is because it is driven by retro-actions between damage localization within a band and the stress field of the surrounding rock that cannot be captured in a pointwise model.

The distinct regimes of crack growth are also responsible for the two stages of creep observed in our constant stress simulations (Figures 4c and 4d, see Section 3), as well as in brittle creep experiments (Figure 1). In the representation of Figure 6a, a constant stress experiment simply maps as a horizontal line starting from any point of the constant strain rate trajectories prior to the peak stress. In order to break the material under constant stress, the $K_I$ trajectory has to remain in the domain $K_I > 0$ up to $D = 1$. There thus exists a threshold in differential stress $\Delta\sigma_{bc}$ that must be met for the sample to fail macroscopically. Otherwise, the accumulation of damage under constant stress will decrease $K_I$ all the way to negative values, inhibiting the growth of further damage. This threshold corresponds to the largest differential stress able to produce a stress concentration factor equal to zero. It is represented in Figure 6a by the summit of the dashed contour of $K_I = 0$, and corresponds to a value of around 720 MPa for a confining pressure of 150 MPa with our inverted Westerly granite parameters. Overall, $\Delta\sigma_{bc}$ can be thought of as a theoretical minimum strength of the rock, which is a function of confining pressure only (continuous blue line in Figure 11). If a brittle creep test is carried out under a constant differential stress above $\Delta\sigma_{bc}$, the damage state will eventually reach $D_c$ in a finite amount of time, then transition to the interacting regime that allows failure. In this case, the creep test will begin by a decrease in $K_I$ that manifests as a decrease in the macroscopic strain rate referred to as decelerating or primary creep (Figure 4c). We note that the minimum brittle
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Brittle yield envelopes of the SCAM model calibrated against Westerly granite data (lines), compared to the experimentally determined strength of intact samples (Byerlee, 1967; Wawersik & Brace, 1971) (squares and circles) and pre-cut samples (Byerlee, 1978) (crosses). Red and orange plain lines correspond to peak strengths at $\dot{e}_{ax} = -10^{-5} \text{ s}^{-1}$ and $-10^{-15} \text{ s}^{-1}$, respectively. The plain blue line corresponds to the minimum failure strength near the limit of infinitely slow strain rates. Modeled failure envelopes exhibit a constant effective friction coefficient of ~0.96 and strain rate-dependent cohesion. The dashed blue line represents the Mohr-Coulomb yield envelope associated with no cohesion and a friction coefficient $\mu = 0.7$, which corresponds to the friction used on shear defects and on macroscopic shear bands in the SCAM simulations.</figcaption></figure>

creep strain rate should be captured accurately in 0-D simulations, since it corresponds to the strain rate at $D \sim D_c$, the extreme value of damage at which damage can grow in a distributed fashion. Finally, as the system switches to the interacting regime, $K_I$ and the macroscopic strain rate both increase, first slowly then catastrophically, accounting for tertiary creep and failure of the sample. The idea that the transition to an accelerating regime of crack interaction up to failure corresponds to crossing a threshold in damage $D_c$ explains the observations of Baud and Meredith (1997), who noted that the transition to tertiary brittle creep coincided with a critical extent of microcracking.

## 5.2. Revisiting the Byerlee Limit

Figure 11 shows failure envelopes for Westerly granite as determined with intact samples (peak stresses at laboratory strain rate, e.g., Byerlee (1967); Wawersik and Brace (1971), black circles and squares), as well as pre-cut samples (the “maximum friction” point from Byerlee (1978), black crosses). Both are on the order of hundreds of MPa, and increase linearly with confining pressures in the 20–200 MPa range. The intact envelope has a steeper slope and lies ~400–500 MPa above the pre-cut strength. Both our 0-D and 2-D models reproduce the intact envelope under the same laboratory strain rate of $10^{-5} \text{ s}^{-1}$ (red line). The simulated envelopes are linear in confining pressure and display an effective cohesion of 94 MPa (inferred by linear regression of the red curve). On the other hand, the standard Mohr-Coulomb plasticity framework, with no cohesion and a friction coefficient of 0.7 provides a good fit to the strength of pre-cut samples (dashed blue line).

If one was to model the transition from intact to broken through strain-softened plasticity (Figure 2a), the friction should drop from ~0.96 to ~0.7, and the cohesion from 94 to 0 MPa. This should occur over a very small amount of plastic strain $\Delta e_{II}^p$ to produce a sharp stress drop. It should be noted that a high “intact” friction coefficient, such as ~0.96, would produce unrealistic shear band orientations (e.g., Coulomb angles of ~23° between the band and $\sigma_1$). Within this model, friction is a property of the bulk material that must evolve as deformation accrues. By contrast, within the SCAM framework, friction is an intrinsic property of planar discontinuities in the rock that manifests at two scales. Friction first conditions slip on small-scale discontinuities (shear defects) whose interaction leads to the formation of larger-scale frictional interfaces (macroscopic shear bands). Those two frictional scales are characterized by the same friction coefficient and no cohesion. Until cracks coalesce, frictional sliding only occurs at the scale of shear defects, and its effect on the material is resolved through its induced stress concentration leading to tensile cracking. After coalescence, the broken material acts as a new frictional zone that generates its own stress perturbations on the surrounding “unbroken” material. It leads, in 2-D setups, to a shear band growing along a direction in which macroscopic stress concentrations amplify damage growth, yielding a
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Brittle yield envelopes of the SCAM model calibrated against Westerly granite data (lines), compared to the experimentally determined strength of intact samples (Byerlee, 1967; Wawersik & Brace, 1971) (squares and circles) and pre-cut samples (Byerlee, 1978) (crosses). Red and orange plain lines correspond to peak strengths at $\dot{e}_{ax} = -10^{-5} \text{ s}^{-1}$ and $-10^{-15} \text{ s}^{-1}$, respectively. The plain blue line corresponds to the minimum failure strength near the limit of infinitely slow strain rates. Modeled failure envelopes exhibit a constant effective friction coefficient of ~0.96 and strain rate-dependent cohesion. The dashed blue line represents the Mohr-Coulomb yield envelope associated with no cohesion and a friction coefficient $\mu = 0.7$, which corresponds to the friction used on shear defects and on macroscopic shear bands in the SCAM simulations.</figcaption></figure>

creep strain rate should be captured accurately in 0-D simulations, since it corresponds to the strain rate at $D \sim D_c$, the extreme value of damage at which damage can grow in a distributed fashion. Finally, as the system switches to the interacting regime, $K_I$ and the macroscopic strain rate both increase, first slowly then catastrophically, accounting for tertiary creep and failure of the sample. The idea that the transition to an accelerating regime of crack interaction up to failure corresponds to crossing a threshold in damage $D_c$ explains the observations of Baud and Meredith (1997), who noted that the transition to tertiary brittle creep coincided with a critical extent of microcracking.

## 5.2. Revisiting the Byerlee Limit

Figure 11 shows failure envelopes for Westerly granite as determined with intact samples (peak stresses at laboratory strain rate, e.g., Byerlee (1967); Wawersik and Brace (1971), black circles and squares), as well as pre-cut samples (the “maximum friction” point from Byerlee (1978), black crosses). Both are on the order of hundreds of MPa, and increase linearly with confining pressures in the 20–200 MPa range. The intact envelope has a steeper slope and lies ~400–500 MPa above the pre-cut strength. Both our 0-D and 2-D models reproduce the intact envelope under the same laboratory strain rate of $10^{-5} \text{ s}^{-1}$ (red line). The simulated envelopes are linear in confining pressure and display an effective cohesion of 94 MPa (inferred by linear regression of the red curve). On the other hand, the standard Mohr-Coulomb plasticity framework, with no cohesion and a friction coefficient of 0.7 provides a good fit to the strength of pre-cut samples (dashed blue line).

If one was to model the transition from intact to broken through strain-softened plasticity (Figure 2a), the friction should drop from ~0.96 to ~0.7, and the cohesion from 94 to 0 MPa. This should occur over a very small amount of plastic strain $\Delta e_{II}^p$ to produce a sharp stress drop. It should be noted that a high “intact” friction coefficient, such as ~0.96, would produce unrealistic shear band orientations (e.g., Coulomb angles of ~23° between the band and $\sigma_1$). Within this model, friction is a property of the bulk material that must evolve as deformation accrues. By contrast, within the SCAM framework, friction is an intrinsic property of planar discontinuities in the rock that manifests at two scales. Friction first conditions slip on small-scale discontinuities (shear defects) whose interaction leads to the formation of larger-scale frictional interfaces (macroscopic shear bands). Those two frictional scales are characterized by the same friction coefficient and no cohesion. Until cracks coalesce, frictional sliding only occurs at the scale of shear defects, and its effect on the material is resolved through its induced stress concentration leading to tensile cracking. After coalescence, the broken material acts as a new frictional zone that generates its own stress perturbations on the surrounding “unbroken” material. It leads, in 2-D setups, to a shear band growing along a direction in which macroscopic stress concentrations amplify damage growth, yielding a
large differential stress drop of hundreds of MPas, which takes place over a very small range of axial strain (Figures 10a and 10b).

It is remarkable that our best fitting coefficient of friction for Westerly granite data (0.7)—which is constrained by data up to the peak strength—also fits the strength of pre-cut samples (dashed blue line in Figure 11). This supports our approach of switching from a damage model to cohesionless Mohr-Coulomb plasticity while retaining a constant coefficient of friction. This approach has the advantage of producing consistent shear band angles of around 30° to the most compressive stress (Figure 10c).

The SCAM framework also allows us to investigate failure at much slower deformation rates because brittle creep data contributes strong constraints on the rate dependence of the pre-peak behavior (Figures 5a and 10b), which is rooted in sub-critical crack growth. As an example, the 0-D failure envelope at a tectonic strain rate of $10^{-15}$ s⁻¹ is shown by the orange line in Figure 11, and represents a constant strength contrast of ~300 MPa relative to laboratory strain rates (red line). Compared to the laboratory strain rate simulations, the tectonic strain rate envelope amounts to lower effective cohesion (~15 MPa) and a similar effective friction. To investigate the model's behavior in the limit of extremely slow strain rates, we construct an estimate of minimum intact strength using the 0-D SCAM model (solid blue line). We calculate it as the differential stress value at $D = D_c$ that would drive a damage growth rate arbitrarily set to 0.1 per billion year. This yields a line with an intercept that is still significantly greater than zero ($C \sim 6.4$ MPa). Overall, the effective cohesion of the brittle failure envelope can be thought of as a strain rate dependent term that does not entirely vanish in the limit of long loading times (e.g., planetary lifetime). The effective friction coefficient however remains invariant with respect to strain rate.

In summary, while standard strain-softened plasticity treats the brittle limit as an envelope that can move with accumulated plastic strain, the SCAM model treats it as a tenuous failure domain (red area in Figure 11) whose upper boundary depends on strain rate. Specifically, the maximum stress that must be attained for the rock to fail macroscopically can vary by ~400 MPa between tectonic and laboratory conditions (10 orders of magnitude in strain rate). This effect is, by definition, not captured by rate-independent elasto-plastic models. On the other hand, the onset of damage growth is not sufficient to define the lower boundary of the failure domain. In order to activate inelastic strain, stresses must exceed a threshold that corresponds to $K_I(D_0) = 0$, that is, the activation of tensile cracking through frictional sliding on small defects. Because frictional sliding is indexed on meeting a cohesionless Mohr-Coulomb criterion with a friction of 0.7, this threshold is closely related to the dashed blue line in Figure 11. Exceeding this threshold however does not guarantee macroscopic failure. If differential stress remains between the dashed and solid blue lines (orange area in Figure 11), damage can be generated in the isolated crack regime, but will never reach $D_c$ in a reasonable amount of time. The rock will thus never fail.

### 5.3. On the Meaning of the SCAM Micromechanical Parameters

Modeling progressive brittle failure with the SCAM model involves a number of parameters that lend themselves to micro-mechanical interpretations. However, the numerous simplifications made by the wing crack model (Ashby & Sammis, 1990), and the use of empirical rules for shear modulus weakening and sub-critical crack growth warrant some caution in doing so. Here we discuss the extent to which parameter values determined by calibrating SCAM against laboratory data provides meaningful information on a rock's microstructure.

Previous studies applying the wing crack micromechanical framework (Ashby & Sammis, 1990), such as Bhat et al. (2011) and Brantut et al. (2012), made the assumption that wing crack initiation happens at $K_I = K_{IC}$. In that case, the analytical expression for $K_I$ (Equation 3 in Ashby and Hallam (1986)) can be used to constrain $\mu$, $a$ and $K_{IC}$ given prior knowledge of $a$ or $K_{IC}$, and of the axial stress of a sample at the onset of microcracking $\sigma_{1c}$, for various confining pressures. The onset of microcracking is typically indexed on the onset of dilatancy or acoustic emissions (Ashby & Sammis, 1990; Brace et al., 1966) (Figure 1). Applying this method to a linear fit of the Brace et al. (1966) data ($\sigma_{1c} = 3.6\sigma_3 + 100$ MPa) yields $\mu = 0.69$, $a = 0.58$ mm, assuming the fracture toughness $K_{IC} = 1.29$ MPa·m$^{1/2}$ determined by our joint inversions in Section 3.3. This coefficient of friction is very close to our estimate, but the crack radius is two orders of magnitude larger than our value of 6.66 µm, closer to the rock's grain size: the length scale considered susceptible to drive wing-cracking. For Darley Dale sandstone, X. Wu et al. (2000) used the same methodology and obtained $\mu = 0.69$ as well as $K_{IC} \le 0.1$ MPa·m$^{1/2}$, assuming $a = 0.11$ mm. These values differ from our inversion results of $\mu = 0.51$, $K_{IC} = 1.4$ MPa·m$^{1/2}$ and $a = 0.65$ mm, especially for $K_{IC}$. These differences may be attributed to our use of a sub-critical crack growth law as opposed to
large differential stress drop of hundreds of MPas, which takes place over a very small range of axial strain (Figures 10a and 10b).

It is remarkable that our best fitting coefficient of friction for Westerly granite data (0.7)—which is constrained by data up to the peak strength—also fits the strength of pre-cut samples (dashed blue line in Figure 11). This supports our approach of switching from a damage model to cohesionless Mohr-Coulomb plasticity while retaining a constant coefficient of friction. This approach has the advantage of producing consistent shear band angles of around 30° to the most compressive stress (Figure 10c).

The SCAM framework also allows us to investigate failure at much slower deformation rates because brittle creep data contributes strong constraints on the rate dependence of the pre-peak behavior (Figures 5a and 10b), which is rooted in sub-critical crack growth. As an example, the 0-D failure envelope at a tectonic strain rate of $10^{-15}$ s⁻¹ is shown by the orange line in Figure 11, and represents a constant strength contrast of ~300 MPa relative to laboratory strain rates (red line). Compared to the laboratory strain rate simulations, the tectonic strain rate envelope amounts to lower effective cohesion (~15 MPa) and a similar effective friction. To investigate the model's behavior in the limit of extremely slow strain rates, we construct an estimate of minimum intact strength using the 0-D SCAM model (solid blue line). We calculate it as the differential stress value at $D = D_c$ that would drive a damage growth rate arbitrarily set to 0.1 per billion year. This yields a line with an intercept that is still significantly greater than zero ($C \sim 6.4$ MPa). Overall, the effective cohesion of the brittle failure envelope can be thought of as a strain rate dependent term that does not entirely vanish in the limit of long loading times (e.g., planetary lifetime). The effective friction coefficient however remains invariant with respect to strain rate.

In summary, while standard strain-softened plasticity treats the brittle limit as an envelope that can move with accumulated plastic strain, the SCAM model treats it as a tenuous failure domain (red area in Figure 11) whose upper boundary depends on strain rate. Specifically, the maximum stress that must be attained for the rock to fail macroscopically can vary by ~400 MPa between tectonic and laboratory conditions (10 orders of magnitude in strain rate). This effect is, by definition, not captured by rate-independent elasto-plastic models. On the other hand, the onset of damage growth is not sufficient to define the lower boundary of the failure domain. In order to activate inelastic strain, stresses must exceed a threshold that corresponds to $K_I(D_0) = 0$, that is, the activation of tensile cracking through frictional sliding on small defects. Because frictional sliding is indexed on meeting a cohesionless Mohr-Coulomb criterion with a friction of 0.7, this threshold is closely related to the dashed blue line in Figure 11. Exceeding this threshold however does not guarantee macroscopic failure. If differential stress remains between the dashed and solid blue lines (orange area in Figure 11), damage can be generated in the isolated crack regime, but will never reach $D_c$ in a reasonable amount of time. The rock will thus never fail.

### 5.3. On the Meaning of the SCAM Micromechanical Parameters

Modeling progressive brittle failure with the SCAM model involves a number of parameters that lend themselves to micro-mechanical interpretations. However, the numerous simplifications made by the wing crack model (Ashby & Sammis, 1990), and the use of empirical rules for shear modulus weakening and sub-critical crack growth warrant some caution in doing so. Here we discuss the extent to which parameter values determined by calibrating SCAM against laboratory data provides meaningful information on a rock's microstructure.

Previous studies applying the wing crack micromechanical framework (Ashby & Sammis, 1990), such as Bhat et al. (2011) and Brantut et al. (2012), made the assumption that wing crack initiation happens at $K_I = K_{IC}$. In that case, the analytical expression for $K_I$ (Equation 3 in Ashby and Hallam (1986)) can be used to constrain $\mu$, $a$ and $K_{IC}$ given prior knowledge of $a$ or $K_{IC}$, and of the axial stress of a sample at the onset of microcracking $\sigma_{1c}$, for various confining pressures. The onset of microcracking is typically indexed on the onset of dilatancy or acoustic emissions (Ashby & Sammis, 1990; Brace et al., 1966) (Figure 1). Applying this method to a linear fit of the Brace et al. (1966) data ($\sigma_{1c} = 3.6\sigma_3 + 100$ MPa) yields $\mu = 0.69$, $a = 0.58$ mm, assuming the fracture toughness $K_{IC} = 1.29$ MPa·m¹/² determined by our joint inversions in Section 3.3. This coefficient of friction is very close to our estimate, but the crack radius is two orders of magnitude larger than our value of 6.66 µm, closer to the rock's grain size: the length scale considered susceptible to drive wing-cracking. For Darley Dale sandstone, X. Wu et al. (2000) used the same methodology and obtained $\mu = 0.69$ as well as $K_{IC} \le 0.1$ MPa·m¹/², assuming $a = 0.11$ mm. These values differ from our inversion results of $\mu = 0.51$, $K_{IC} = 1.4$ MPa·m¹/² and $a = 0.65$ mm, especially for $K_{IC}$. These differences may be attributed to our use of a sub-critical crack growth law as opposed to