$K_I = K_{IC}$. We also note that because $K_I$ scales as $\sqrt{a}$, our inversions do a better job at constraining a ratio of $K_{IC}/\sqrt{a}$ rather than $K_I$ or $a$ individually. Some trade-off is therefore expected between these values. It is noteworthy that the data against which the SCAM model parameters were inverted are much more complete than the data used in the studies discussed above. Specifically, our data includes stress-strain curves up to the peak stress (as opposed to a single point: the onset of dilatancy). We however acknowledge that at the peak stress, the discrepancy between the wing-crack model and the real rock microstructure may become large, introducing some bias in our parameter estimate.

Inspecting Figure 12 of Ashby and Sammis (1990) also reveals that wing-crack based damage mechanics can fit the very first portion of the experimental failure envelope of Westerly granite (i.e., at confining pressures between 0 and a few 100 s of MPas), but fail to predict peak stresses at greater pressures due to the non-linear dependence of intact rock strength with respect to confining pressure (visible in Figure 11: squares and circles). Ashby and Sammis (1990) suggested that the curvature of the experimental failure envelope of Westerly granite might result from low temperature ductile flow occurring within weak granite minerals. This hypothesis was tested by Bhat et al. (2011) by modeling a bi-mineralic quartz-feldspar assemblage with a dislocation glide flow law. It was found consistent with the experimental failure envelope at greater confining pressures. Broadly speaking, inverting experimental data over a large range of pressures using a model based only on fracture growth likely neglects processes that can significantly affect the mechanical response, such as intra-grain dislocation glide. This leads to inverted fracture parameters whose value can deviate from their expected range, since they are forced to explain behaviors caused by other deformation mechanisms. The exact value of the parameters should therefore not be over-interpreted. The value of the SCAM model lies more in its ability to extrapolate micromechanics-based rock behavior to larger scales, rather than in its informative power about rocks' intrinsic parameters.

# 6. Toward an Application to Tectonic Problems

In Section 3, we introduced and calibrated the 0-D SCAM micro-mechanical model against experimental data to capture the mechanical behaviors of Westerly granite and Darley Dale sandstone in the brittle regime. Section 4 was dedicated to the implementation of this model in a 2-D plane strain tectonic solver to investigate the behavior of the calibrated model when strain localization is made possible. We validated the model by showing that it still accurately predicts experimental data in 2-D. Next, we turn to the end-goal of the SCAM framework by showcasing initial attempts at using it to model tectonic deformation.

## 6.1. Crustal-Scale Numerical Setup

As a first step toward using the SCAM model for long-term tectonic problems, we focus on the initial stages of faulting of a 10-km thick brittle plate overlying a 40-km wide and 10-km thick low-viscosity Newtonian medium subjected to gravity (9.8 m s⁻²) and to a constant horizontal extension rate (Figure 12). An additional upper layer of low-viscosity “sticky air” ensures the brittle plate has a traction-free top boundary (Gerya, 2010). The plate is notched at the center of its lower edge by a 1 km long and 0.75 km thick protrusion of the underlying viscous layer to promote strain localization in the middle of the domain. The latter is discretized using a cell size of 100 × 250 m within 8 km of the top and bottom walls, and an greater resolution, with a cell size of 100 × 100 m in the remainder of the domain containing the brittle plate. The plate and the fluid layer underneath it are assigned a density of 2,700 kg m⁻³, while the uppermost fluid layer's density is 0.01 kg m⁻³, to ensure negligible pressures at the top of the plate. The left and right sides are prescribed a fixed outward horizontal velocity amounting to a constant extension rate of 10 cm/yr (strain rate of ~10⁻¹³ s⁻¹). The velocities of the top and bottom boundaries are set to satisfy volume conservation within the domain and to preserve the height of the brittle plate's top surface. All boundaries are free slip. As in simulated experiments, the maximum viscosity ηmax allowed in the domain is chosen such that its Maxwell time (η/G) is 50 times longer than the longest simulation time. This guarantees that the brittle plate—whose viscosity is set to ηmax—retains an elastic response throughout the entire simulation. The exact value of ηmax will therefore not matter as long as it is large enough. To avoid large, computationally challenging viscosity contrasts within the model domain, the lower bound on viscosity is set to six orders of magnitude below the brittle plate's viscosity, and is assigned to the air and viscous lower layer. This guarantees that they behave as low-viscosity fluids during each simulation. This setup resembles that used by Lavier et al. (2000) and Olive et al. (2016) to investigate the effect of brittle strain softening and elasticity on extensional tectonic styles, within the standard elasto-plastic framework. One notable difference is that here the base of the
$K_I = K_{IC}$. We also note that because $K_I$ scales as $\sqrt{a}$, our inversions do a better job at constraining a ratio of $K_{IC}/\sqrt{a}$ rather than $K_I$ or $a$ individually. Some trade-off is therefore expected between these values. It is noteworthy that the data against which the SCAM model parameters were inverted are much more complete than the data used in the studies discussed above. Specifically, our data includes stress-strain curves up to the peak stress (as opposed to a single point: the onset of dilatancy). We however acknowledge that at the peak stress, the discrepancy between the wing-crack model and the real rock microstructure may become large, introducing some bias in our parameter estimate.

Inspecting Figure 12 of Ashby and Sammis (1990) also reveals that wing-crack based damage mechanics can fit the very first portion of the experimental failure envelope of Westerly granite (i.e., at confining pressures between 0 and a few 100 s of MPas), but fail to predict peak stresses at greater pressures due to the non-linear dependence of intact rock strength with respect to confining pressure (visible in Figure 11: squares and circles). Ashby and Sammis (1990) suggested that the curvature of the experimental failure envelope of Westerly granite might result from low temperature ductile flow occurring within weak granite minerals. This hypothesis was tested by Bhat et al. (2011) by modeling a bi-mineralic quartz-feldspar assemblage with a dislocation glide flow law. It was found consistent with the experimental failure envelope at greater confining pressures. Broadly speaking, inverting experimental data over a large range of pressures using a model based only on fracture growth likely neglects processes that can significantly affect the mechanical response, such as intra-grain dislocation glide. This leads to inverted fracture parameters whose value can deviate from their expected range, since they are forced to explain behaviors caused by other deformation mechanisms. The exact value of the parameters should therefore not be over-interpreted. The value of the SCAM model lies more in its ability to extrapolate micromechanics-based rock behavior to larger scales, rather than in its informative power about rocks' intrinsic parameters.

# 6. Toward an Application to Tectonic Problems

In Section 3, we introduced and calibrated the 0-D SCAM micro-mechanical model against experimental data to capture the mechanical behaviors of Westerly granite and Darley Dale sandstone in the brittle regime. Section 4 was dedicated to the implementation of this model in a 2-D plane strain tectonic solver to investigate the behavior of the calibrated model when strain localization is made possible. We validated the model by showing that it still accurately predicts experimental data in 2-D. Next, we turn to the end-goal of the SCAM framework by showcasing initial attempts at using it to model tectonic deformation.

## 6.1. Crustal-Scale Numerical Setup

As a first step toward using the SCAM model for long-term tectonic problems, we focus on the initial stages of faulting of a 10-km thick brittle plate overlying a 40-km wide and 10-km thick low-viscosity Newtonian medium subjected to gravity (9.8 m s⁻²) and to a constant horizontal extension rate (Figure 12). An additional upper layer of low-viscosity “sticky air” ensures the brittle plate has a traction-free top boundary (Gerya, 2010). The plate is notched at the center of its lower edge by a 1 km long and 0.75 km thick protrusion of the underlying viscous layer to promote strain localization in the middle of the domain. The latter is discretized using a cell size of 100 × 250 m within 8 km of the top and bottom walls, and an greater resolution, with a cell size of 100 × 100 m in the remainder of the domain containing the brittle plate. The plate and the fluid layer underneath it are assigned a density of 2,700 kg m⁻³, while the uppermost fluid layer's density is 0.01 kg m⁻³, to ensure negligible pressures at the top of the plate. The left and right sides are prescribed a fixed outward horizontal velocity amounting to a constant extension rate of 10 cm/yr (strain rate of ~10⁻¹³ s⁻¹). The velocities of the top and bottom boundaries are set to satisfy volume conservation within the domain and to preserve the height of the brittle plate's top surface. All boundaries are free slip. As in simulated experiments, the maximum viscosity ηmax allowed in the domain is chosen such that its Maxwell time (η/G) is 50 times longer than the longest simulation time. This guarantees that the brittle plate—whose viscosity is set to ηmax—retains an elastic response throughout the entire simulation. The exact value of ηmax will therefore not matter as long as it is large enough. To avoid large, computationally challenging viscosity contrasts within the model domain, the lower bound on viscosity is set to six orders of magnitude below the brittle plate's viscosity, and is assigned to the air and viscous lower layer. This guarantees that they behave as low-viscosity fluids during each simulation. This setup resembles that used by Lavier et al. (2000) and Olive et al. (2016) to investigate the effect of brittle strain softening and elasticity on extensional tectonic styles, within the standard elasto-plastic framework. One notable difference is that here the base of the
<figure><img src="image_14.png" /></figure>

**Figure 12.** Numerical setup used to simulate the stretching of a rigid crustal unit in 2-D plane strain. The domain is 40 km long and 30 km thick and contains two 10-km thick low-viscosity layers at the top and bottom of a brittle plate of the same thickness. The brittle plate is notched by a 1 × 0.75 km protrusion of the underlying viscous fluid. Constant outward velocities at the left and right boundaries apply a constant stretching rate on the brittle plate, while the top and bottom boundaries are also assigned constant inward velocities to satisfy volume conservation. All boundaries are free of shear tractions.

brittle plate is a lithological boundary that gets advected as the plate thins (unlike, e.g., a thermal boundary that may experience diffusion).

We run two suites of simulations, one with the SCAM parameterization of brittle failure in Westerly granite (Table 1), and one with the standard strain weakened elasto-plastic (EP) approach, similar to the runs of Olive et al. (2016). Tectonic simulations using the SCAM model, like those presented in Section 4, include a switch to standard elasto-plasticity in fully damaged areas. Specifically, the damage viscosity is set to smoothly transition to a plastic viscosity as $\eta_D$ approaches $\eta_p$ over a viscosity range of $|\eta_D - \eta_p|$ equal to $\eta_{min}/50$. The locally broken material (i.e., $D = 1$) behaves as a Mohr-Coulomb plastic solid, with no cohesion and the same friction coefficient as that of the shear defects.

Standard elasto-plastic simulations, on the other hand, are parameterized to match the intact strength of the SCAM model calibrated on Westerly granite (Figure 11) under a laboratory strain rate ($10^{-5} s^{-1}$, red curve) and under a tectonic strain rate ($10^{-15} s^{-1}$, orange curve). These envelopes correspond to a similar friction coefficient $\mu = 0.96$ and cohesions $C$ of 94 and 15 MPa, respectively (see Section 5.2). For concision, these simulations will be referred to as EP-lab and EP-tecto. In both cases the frictional parameters are linearly weakened over a critical amount of accumulated plastic strain $\Delta\epsilon_{II} = 0.1$ down to $\mu = 0.7$ and $C = 0$ MPa, the frictional properties of the shear defects.

## 6.2. Development of Fault Networks

We first compare the faulting patterns produced by an elasto-plastic rheology (EP-tecto) versus the SCAM rheology during early rifting (up to ~180 m of total extension). Figure 13 shows successive snapshots of accumulated plastic strain $\epsilon_{II}^p$ (a1 to a3), second invariant of strain rate $\dot{e}_{II}$ (b1 to b3), and viscosity (c1 to c3), as extension of the EP-tecto plate progresses. Plastic yielding starts from the surface of the plate, where pressure, and therefore yield stress, is lowest. It progressively deepens as the yield criterion is met deeper and deeper due to elastic loading of the plate (Panels a1 and a2). Plastic yielding initially mostly develops in a distributed fashion
<figure><img src="image_14.png" /></figure>

**Figure 12.** Numerical setup used to simulate the stretching of a rigid crustal unit in 2-D plane strain. The domain is 40 km long and 30 km thick and contains two 10-km thick low-viscosity layers at the top and bottom of a brittle plate of the same thickness. The brittle plate is notched by a 1 × 0.75 km protrusion of the underlying viscous fluid. Constant outward velocities at the left and right boundaries apply a constant stretching rate on the brittle plate, while the top and bottom boundaries are also assigned constant inward velocities to satisfy volume conservation. All boundaries are free of shear tractions.

brittle plate is a lithological boundary that gets advected as the plate thins (unlike, e.g., a thermal boundary that may experience diffusion).

We run two suites of simulations, one with the SCAM parameterization of brittle failure in Westerly granite (Table 1), and one with the standard strain weakened elasto-plastic (EP) approach, similar to the runs of Olive et al. (2016). Tectonic simulations using the SCAM model, like those presented in Section 4, include a switch to standard elasto-plasticity in fully damaged areas. Specifically, the damage viscosity is set to smoothly transition to a plastic viscosity as $\eta_D$ approaches $\eta_p$ over a viscosity range of $|\eta_D - \eta_p|$ equal to $\eta_{min}/50$. The locally broken material (i.e., $D=1$) behaves as a Mohr-Coulomb plastic solid, with no cohesion and the same friction coefficient as that of the shear defects.

Standard elasto-plastic simulations, on the other hand, are parameterized to match the intact strength of the SCAM model calibrated on Westerly granite (Figure 11) under a laboratory strain rate ($10^{-5} s^{-1}$, red curve) and under a tectonic strain rate ($10^{-15} s^{-1}$, orange curve). These envelopes correspond to a similar friction coefficient $\mu = 0.96$ and cohesions $C$ of 94 and 15 MPa, respectively (see Section 5.2). For concision, these simulations will be referred to as EP-lab and EP-tecto. In both cases the frictional parameters are linearly weakened over a critical amount of accumulated plastic strain $\Delta\epsilon_{II} = 0.1$ down to $\mu = 0.7$ and $C = 0$ MPa, the frictional properties of the shear defects.

## 6.2. Development of Fault Networks

We first compare the faulting patterns produced by an elasto-plastic rheology (EP-tecto) versus the SCAM rheology during early rifting (up to ~180 m of total extension). Figure 13 shows successive snapshots of accumulated plastic strain $\epsilon_{II}^p$ (a1 to a3), second invariant of strain rate $\dot{e}_{II}$ (b1 to b3), and viscosity (c1 to c3), as extension of the EP-tecto plate progresses. Plastic yielding starts from the surface of the plate, where pressure, and therefore yield stress, is lowest. It progressively deepens as the yield criterion is met deeper and deeper due to elastic loading of the plate (Panels a1 and a2). Plastic yielding initially mostly develops in a distributed fashion
<figure><img src="image_15.png" /><figcaption>**Figure 13.** Snapshots of an early rifting simulation (10 cm/yr) in a 10-km thick brittle plate using the elasto-plastic EP-tecto parameterization (a–c) versus the SCAM model (d–f). Each row corresponds to a specific time (i.e., amount of finite extension) in the simulation. The first column shows the accumulated plastic strain for EP-tecto, and the damage field for the SCAM model. The second column shows the second invariant of the deviatoric strain rate tensor, and the third column displays the viscosity field within the brittle plate.</figcaption></figure>

(second row), with plastic strain localization occurring once almost half of the plate has reached yielding. The zone of distributed yielding closer to the viscous protrusion, then transforms into an area densely populated with shear bands, of dip angle near 55°, each accommodating a very small fraction of the total extension rate (second row). Stress concentrations around the basal notch eventually lead to the formation of a pair of shear bands symmetrically cutting across the brittle plate, following the path of pre-existing superficial shear bands. Plastic strain accordingly localizes along two major antithetic shear bands (third row). The formation of these faults relax the elastic stresses within the plate and inhibit plastic yielding in the remainder of the plate (Panels b3 and c3).

The same setup using the SCAM flow law (Figures 13d–13f) shows a different story. Damage first increases uniformly within ~1 km below the surface of the brittle plate (Panel d1, blue arrow 1). After shallow damage exceeds values of 0.4–0.5 (~D_c), damage localization proceeds through the downward propagation of damage bands which promptly turn into plastic shear bands (first row). This occurs when the isolated crack regime transitions to the interacting cracks regime (see Section 5.1.2). The thickness of the shallow distributed damage zone (~1 km) is set by a competition between damage growth being activated deeper and deeper as K_I turns positive with loading, and the unstable growth of damage bands (which soon turns plastic) once D exceeds D_c. Damage banding unloads the surrounding material by releasing elastic stress through inelastic flow. The growth of these shear bands is self-promoted as the result of stress concentrations at their lower tips, which locally accelerate crack growth and coalescence. This dynamic is clearly visible on the viscosity snapshots, for example, Panel F1, where each fault has its own associated lobe of low viscosity (blue arrow 2), driven by damage growth. These lobes are accompanied by a front of low (damage) viscosity which deepens through time (Panel f2, blue arrow 3). This reflects the fact that with continued far-field loading, the depth at which K_I turns positive increases.

Growing faults tend to shield each other and alter the stress field in their vicinity, which leads to the development of complex networks shaped in a tree-like fashion (e.g., the imbricated fault structure that forms on the left side of
<figure><img src="image_15.png" /><figcaption>**Figure 13.** Snapshots of an early rifting simulation (10 cm/yr) in a 10-km thick brittle plate using the elasto-plastic EP-tecto parameterization (a–c) versus the SCAM model (d–f). Each row corresponds to a specific time (i.e., amount of finite extension) in the simulation. The first column shows the accumulated plastic strain for EP-tecto, and the damage field for the SCAM model. The second column shows the second invariant of the deviatoric strain rate tensor, and the third column displays the viscosity field within the brittle plate.</figcaption></figure>

(second row), with plastic strain localization occurring once almost half of the plate has reached yielding. The zone of distributed yielding closer to the viscous protrusion, then transforms into an area densely populated with shear bands, of dip angle near 55°, each accommodating a very small fraction of the total extension rate (second row). Stress concentrations around the basal notch eventually lead to the formation of a pair of shear bands symmetrically cutting across the brittle plate, following the path of pre-existing superficial shear bands. Plastic strain accordingly localizes along two major antithetic shear bands (third row). The formation of these faults relax the elastic stresses within the plate and inhibit plastic yielding in the remainder of the plate (Panels b3 and c3).

The same setup using the SCAM flow law (Figures 13d–13f) shows a different story. Damage first increases uniformly within ~1 km below the surface of the brittle plate (Panel d1, blue arrow 1). After shallow damage exceeds values of 0.4–0.5 (~D_c), damage localization proceeds through the downward propagation of damage bands which promptly turn into plastic shear bands (first row). This occurs when the isolated crack regime transitions to the interacting cracks regime (see Section 5.1.2). The thickness of the shallow distributed damage zone (~1 km) is set by a competition between damage growth being activated deeper and deeper as K_I turns positive with loading, and the unstable growth of damage bands (which soon turns plastic) once D exceeds D_c. Damage banding unloads the surrounding material by releasing elastic stress through inelastic flow. The growth of these shear bands is self-promoted as the result of stress concentrations at their lower tips, which locally accelerate crack growth and coalescence. This dynamic is clearly visible on the viscosity snapshots, for example, Panel F1, where each fault has its own associated lobe of low viscosity (blue arrow 2), driven by damage growth. These lobes are accompanied by a front of low (damage) viscosity which deepens through time (Panel f2, blue arrow 3). This reflects the fact that with continued far-field loading, the depth at which K_I turns positive increases.

Growing faults tend to shield each other and alter the stress field in their vicinity, which leads to the development of complex networks shaped in a tree-like fashion (e.g., the imbricated fault structure that forms on the left side of
the plate in Panels d3–f3). After about 50 m of stretching, the concentration of tensile stresses around the viscous protrusion promotes the growth of a major fault, which first grows upward with a downward convex shape (Figures 13d2–13f2), connecting with one of the deepening faults (Panel d2, blue arrow 4), before a more favorably oriented branch eventually grows and bypasses the less favorably oriented pre-existing segment of the fault. This master-fault's dip angle increases from ~60° at the bottom of the plate, to almost 90° close to the surface. As extension progresses, the uppermost part of the fault (Panel d3, blue arrow 5) rotates toward gentler dips, which progressively widens the thickness of the damage band near the surface, as shown in the last row of snapshots (blue arrow 5 in Panel D3). The dip angle of secondary faults is variable and ranges from 60° to 75° with segments locally close to vertical, reflecting the spatially variable distribution of stresses induced by the complex array of faults.

The fault angles with respect to $σ_1$ (here vertical) are overall smaller in the SCAM tectonic simulations (15–30°) than what was observed in our 2-D simulations of triaxial experiments (~30°, see Section 4.4). These faults are closer to the Coulomb angle (23°) associated with the effective intact strength coefficient of internal friction of Westerly Granite ($μ$ = 0.96), than with the coefficient of friction at the micro-crack scale ($μ$ = 0.7). Some deviation from the theoretical value likely stems from heterogeneities in the stress field that develop through stress concentrations and complex interactions between growing shear bands. In simulation EP-tecto, the master fault orientations of around 33°, correspond to the Arthur angle $θ_A$ = 45° – (φ + $ψ$)/4 (Arthur & Dunstan, 1977) associated with the prescribed initial internal friction coefficient of 0.96 without any plastic dilatancy ($ψ$ = 0), which is a typical shear band angle found (along with the Coulomb angle) in numerical elasto-plastic simulations (Kaus, 2010).

Contrary to the EP-tecto simulation, the initiation of a fault cutting across the plate does not inhibit the activity on secondary faults, which continue to accommodate a significant amount of extension in the SCAM run (Figure 13, Panel b3 vs. e3). We note that EP simulations can produce such behavior, but it generally requires a small strength contrast between the fault zone and surrounding lithosphere (Lavier et al., 2000). To further characterize the partitioning of strain in the SCAM runs, Figure 14a represents the distribution of slip rates on the population of faults, sampled at different depths in the plate after ~180 m of extension (Figures 13a3–13c3 and 13d3–13f3). We identify normal faults by locating positive peaks along a profile of horizontal strain rate $\dot{e}_{xx}$ along a line of constant depth. To mitigate resolution issues when two neighboring peaks are found, we only retain the larger of any two peaks distant by less than four cell sizes. We focus on peaks where strain rate exceeds $4 \cdot 10^{-13} s^{-1}$, which

<figure><img src="image_16.png" /><figcaption>**Figure 14.** (a) Distribution of slip rates on fault populations sampled along different depths in a SCAM-based simulation (gray lines) versus simulation EP-tecto (yellow-to-red lines). Faults are measured after 180 m of total stretching, corresponding to Panels a3–c3 and d3–f3 of Figure 13. (b) Differential stress averaged across a 1 km-wide vertical band along the left side of the extended brittle plate, plotted as a function of total extensional strain. The black curve corresponds to the SCAM simulation, whereas the red and orange curves are associated to elasto-plastic simulations parameterized with short (EP-lab) and long-term (EP-tecto) intact strengths. Stars mark the timing of the snapshots displayed in Figure 13.</figcaption></figure>
the plate in Panels d3–f3). After about 50 m of stretching, the concentration of tensile stresses around the viscous protrusion promotes the growth of a major fault, which first grows upward with a downward convex shape (Figures 13d2–13f2), connecting with one of the deepening faults (Panel d2, blue arrow 4), before a more favorably oriented branch eventually grows and bypasses the less favorably oriented pre-existing segment of the fault. This master-fault's dip angle increases from ~60° at the bottom of the plate, to almost 90° close to the surface. As extension progresses, the uppermost part of the fault (Panel d3, blue arrow 5) rotates toward gentler dips, which progressively widens the thickness of the damage band near the surface, as shown in the last row of snapshots (blue arrow 5 in Panel D3). The dip angle of secondary faults is variable and ranges from 60° to 75° with segments locally close to vertical, reflecting the spatially variable distribution of stresses induced by the complex array of faults.

The fault angles with respect to $σ_1$ (here vertical) are overall smaller in the SCAM tectonic simulations (15–30°) than what was observed in our 2-D simulations of triaxial experiments (~30°, see Section 4.4). These faults are closer to the Coulomb angle (23°) associated with the effective intact strength coefficient of internal friction of Westerly Granite ($μ$ = 0.96), than with the coefficient of friction at the micro-crack scale ($μ$ = 0.7). Some deviation from the theoretical value likely stems from heterogeneities in the stress field that develop through stress concentrations and complex interactions between growing shear bands. In simulation EP-tecto, the master fault orientations of around 33°, correspond to the Arthur angle $θ_A$ = 45° – (φ + $ψ$)/4 (Arthur & Dunstan, 1977) associated with the prescribed initial internal friction coefficient of 0.96 without any plastic dilatancy ($ψ$ = 0), which is a typical shear band angle found (along with the Coulomb angle) in numerical elasto-plastic simulations (Kaus, 2010).

Contrary to the EP-tecto simulation, the initiation of a fault cutting across the plate does not inhibit the activity on secondary faults, which continue to accommodate a significant amount of extension in the SCAM run (Figure 13, Panel b3 vs. e3). We note that EP simulations can produce such behavior, but it generally requires a small strength contrast between the fault zone and surrounding lithosphere (Lavier et al., 2000). To further characterize the partitioning of strain in the SCAM runs, Figure 14a represents the distribution of slip rates on the population of faults, sampled at different depths in the plate after ~180 m of extension (Figures 13a3–13c3 and 13d3–13f3). We identify normal faults by locating positive peaks along a profile of horizontal strain rate $\dot{e}_{xx}$ along a line of constant depth. To mitigate resolution issues when two neighboring peaks are found, we only retain the larger of any two peaks distant by less than four cell sizes. We focus on peaks where strain rate exceeds $4 \cdot 10^{-13} s^{-1}$, which

<figure><img src="image_16.png" /><figcaption>**Figure 14.** (a) Distribution of slip rates on fault populations sampled along different depths in a SCAM-based simulation (gray lines) versus simulation EP-tecto (yellow-to-red lines). Faults are measured after 180 m of total stretching, corresponding to Panels a3–c3 and d3–f3 of Figure 13. (b) Differential stress averaged across a 1 km-wide vertical band along the left side of the extended brittle plate, plotted as a function of total extensional strain. The black curve corresponds to the SCAM simulation, whereas the red and orange curves are associated to elasto-plastic simulations parameterized with short (EP-lab) and long-term (EP-tecto) intact strengths. Stars mark the timing of the snapshots displayed in Figure 13.</figcaption></figure>
amounts to five times the background horizontal strain rate imposed by boundary conditions. The slip rate on each fault is then calculated assuming that the thickness of the faults is one cell length (100 m in these simulations).

Lastly, we find that the slip rates of the incipient fault population produced by the SCAM model follows a power law distribution of exponent close to 1 (gray lines in Figure 14a). This distribution likely originates in the tree-like, near-fractal nature of the fault network that emerges within the brittle plate. Simulation EP-tecto, on the other hand, rapidly produces two dominant faults with very similar slip rates, and virtually no smaller fault with slower slip rates. The distribution of fault slip rates produced by the SCAM run evokes the power law distributions exhibited by natural populations of normal faults. For example, the lengths and offsets of Basin and Range normal faults have been shown to follow a power law distribution with exponent close to unity (e.g., Marrett et al., 1999; Scholz et al., 1993). Natural rift systems are also known to partition extension onto multiple faults of varying sizes, such that all minor intrabasinal faults typically accommodate as much strain as a single half-graben border fault (Morley, 1996). Standard elasto-plastic models tend to localize deformation onto a small number of lithosphere-scale faults (Olive et al., 2016), but can also produce a wide range of fault sizes in the early stages of rifting, provided sufficient resolution is used (Naliboff et al., 2020; Pan et al., 2022). In this case, the distribution of fault sizes and offsets appears to depend on the amount and rate of strain softening, as well as on the overall extension rate (Naliboff et al., 2020). The SCAM model provides a novel framework to further explore this dependence without resorting to ad-hoc softening rules.

## 6.3. Strength of the Brittle Lithosphere

Figure 14b shows the differential stress averaged along the left side of the model domain throughout our early rifting simulations. All show an initial elastic increase, followed by a peak and a softening phase. The average peak stress reached in the SCAM simulation is ~100 MPa. The softening phase is much less acute than in our simulations of triaxial experiments (Figure 10). This possibly reflects the progressive nature of the development of the fault network (Figure 13). The peak stress roughly coincides with the moment when a fault first connects the top to the bottom of the plate, shortly after snapshots D2–F2 in Figure 13. The steady, long-term stress (~90 MPa in Figure 14b) is attained when one of the throughcutting faults has clearly developed into the dominant one accommodating the largest extension rate (Figure 13e3).

On the other hand, simulation EP-lab illustrates what happens when one models rifting with strain-softened elasto-plasticity calibrated on experiments conducted at laboratory strain rates (Figure 11). The plate reaches a significantly greater peak stress of ~160 MPa (Figure 14b), and experiences more drastic weakening down to its long-term frictional strength, which is inferred from pre-cut samples ($\mu$ = 0.7 and no cohesion). By contrast, simulation EP-tecto, which is calibrated to reproduce the intact strength of samples predicted by the SCAM model at tectonic strain rates (yellow line in Figure 11), unsurprisingly produces a stress-strain curve similar to the SCAM simulation (Figure 14b).

Through its built-in dependence on strain rate, the SCAM rheology can extrapolate crustal strength under tectonic conditions, even though it is entirely calibrated on laboratory data. It thus constitutes a promising alternative to the standard approach in tectonic modeling, which consists of using “Byerlee’s law,” that is, assigning the pre-cut sample strength to intact brittle lithosphere ((Brace & Kohlstedt, 1980)). The SCAM flow law has the advantage of using a single value of friction and a handful of micro-mechanical parameters to predict the lithosphere’s intact strength across 10 orders of magnitude of strain rate. It self-consistently handles the transition from intact to broken (pre-cut), and does not require any empirical assumption on weakening strain. Alternatively, SCAM can be used to prescribe an intact strength that is appropriate at tectonic strain rates, for use in standard elasto-plastic models (e.g., simulation EP-tecto). This approach, however, does not solve the problem of the ad-hoc weakening strain.

## 6.4. Growth of Geological Structures

In order to compare the tectonic structures formed after greater amounts of finite extension by the SCAM and elasto-plastic (EP-tecto) rheologies, we performed three additional simulations under strain rates of 1, 1, and 10 mm/yr for each rheology. Due to the computationally demanding nature of the SCAM model, we decreased the resolution of the spatial domain using a cell size of 250 × 500 m within 8 km of the top and bottom walls, and a greater resolution, with a cell size of 250 × 250 m in the rest of the domain containing the brittle plate.
amounts to five times the background horizontal strain rate imposed by boundary conditions. The slip rate on each fault is then calculated assuming that the thickness of the faults is one cell length (100 m in these simulations).

Lastly, we find that the slip rates of the incipient fault population produced by the SCAM model follows a power law distribution of exponent close to 1 (gray lines in Figure 14a). This distribution likely originates in the tree-like, near-fractal nature of the fault network that emerges within the brittle plate. Simulation EP-tecto, on the other hand, rapidly produces two dominant faults with very similar slip rates, and virtually no smaller fault with slower slip rates. The distribution of fault slip rates produced by the SCAM run evokes the power law distributions exhibited by natural populations of normal faults. For example, the lengths and offsets of Basin and Range normal faults have been shown to follow a power law distribution with exponent close to unity (e.g., Marrett et al., 1999; Scholz et al., 1993). Natural rift systems are also known to partition extension onto multiple faults of varying sizes, such that all minor intrabasinal faults typically accommodate as much strain as a single half-graben border fault (Morley, 1996). Standard elasto-plastic models tend to localize deformation onto a small number of lithosphere-scale faults (Olive et al., 2016), but can also produce a wide range of fault sizes in the early stages of rifting, provided sufficient resolution is used (Naliboff et al., 2020; Pan et al., 2022). In this case, the distribution of fault sizes and offsets appears to depend on the amount and rate of strain softening, as well as on the overall extension rate (Naliboff et al., 2020). The SCAM model provides a novel framework to further explore this dependence without resorting to ad-hoc softening rules.

## 6.3. Strength of the Brittle Lithosphere

Figure 14b shows the differential stress averaged along the left side of the model domain throughout our early rifting simulations. All show an initial elastic increase, followed by a peak and a softening phase. The average peak stress reached in the SCAM simulation is ~100 MPa. The softening phase is much less acute than in our simulations of triaxial experiments (Figure 10). This possibly reflects the progressive nature of the development of the fault network (Figure 13). The peak stress roughly coincides with the moment when a fault first connects the top to the bottom of the plate, shortly after snapshots D2–F2 in Figure 13. The steady, long-term stress (~90 MPa in Figure 14b) is attained when one of the throughcutting faults has clearly developed into the dominant one accommodating the largest extension rate (Figure 13e3).

On the other hand, simulation EP-lab illustrates what happens when one models rifting with strain-softened elasto-plasticity calibrated on experiments conducted at laboratory strain rates (Figure 11). The plate reaches a significantly greater peak stress of ~160 MPa (Figure 14b), and experiences more drastic weakening down to its long-term frictional strength, which is inferred from pre-cut samples ($\mu$ = 0.7 and no cohesion). By contrast, simulation EP-tecto, which is calibrated to reproduce the intact strength of samples predicted by the SCAM model at tectonic strain rates (yellow line in Figure 11), unsurprisingly produces a stress-strain curve similar to the SCAM simulation (Figure 14b).

Through its built-in dependence on strain rate, the SCAM rheology can extrapolate crustal strength under tectonic conditions, even though it is entirely calibrated on laboratory data. It thus constitutes a promising alternative to the standard approach in tectonic modeling, which consists of using “Byerlee’s law,” that is, assigning the pre-cut sample strength to intact brittle lithosphere ((Brace & Kohlstedt, 1980)). The SCAM flow law has the advantage of using a single value of friction and a handful of micro-mechanical parameters to predict the lithosphere’s intact strength across 10 orders of magnitude of strain rate. It self-consistently handles the transition from intact to broken (pre-cut), and does not require any empirical assumption on weakening strain. Alternatively, SCAM can be used to prescribe an intact strength that is appropriate at tectonic strain rates, for use in standard elasto-plastic models (e.g., simulation EP-tecto). This approach, however, does not solve the problem of the ad-hoc weakening strain.

## 6.4. Growth of Geological Structures

In order to compare the tectonic structures formed after greater amounts of finite extension by the SCAM and elasto-plastic (EP-tecto) rheologies, we performed three additional simulations under strain rates of 1, 1, and 10 mm/yr for each rheology. Due to the computationally demanding nature of the SCAM model, we decreased the resolution of the spatial domain using a cell size of 250 × 500 m within 8 km of the top and bottom walls, and a greater resolution, with a cell size of 250 × 250 m in the rest of the domain containing the brittle plate.
Snapshots of the second invariant of strain rate ($\dot{\epsilon}_{II}$) fields at different times are displayed in Figure 15. The first three rows correspond to successive times throughout the SCAM simulations after ~251, ~1,079, and ~3,457 m of horizontal stretching, respectively. The fourth row corresponds to snapshots of elasto-plastic simulations performed under the same extension rates, after ~3,457 m of extension. The elasto-plastic simulations remain unsurprinsingly insensitive to a change of strain rate (e.g., Olive et al., 2016), as they all develop the same central graben structure (Figures 15d–15f).

The SCAM simulations, on the other hand, generate features that vary with extension rate. The first one is the location of the cluster of faults in the hanging-wall (Figure 15a1–15c1). The faster the extension rate, the closer this cluster of secondary faults lies to the master fault, with distances ranging from 15 to 7 km. The second feature is the geometry of the rift after 3,457 m of extension (Panels a3–c3). The rift structure shows a greater degree of asymmetry with increasing extension rate. This is particularly well expressed in topography, which corresponds to a symmetric graben in Panel a3 but is closer to a half-graben in Panels b3 and c3, as more extension is accommodated of the central right-dipping fault. We acknowledge that some of this variability could be attributed to stochasticity in marker positions, coupled with strong non-linearities of the SCAM flow law. We however ran the SCAM simulations a second time using different randomly assigned marker positions, and otherwise identical parameters. This second set yielded the same pattern of strain rate dependence as illustrated in Panels a3–c3.

Extension rate is known to influence rifting styles, primarily by modulating the thermal structure (Buck, 1991; Lavier & Buck, 2002). Our results suggest that it may play an additional, more subtle role by modulating the very processes of fault development. Specifically, the SCAM model introduces a strain rate dependence of brittle deformation by transiently activating a moderate to low viscosity in portions of the upper crust where damage is actively growing (e.g., Panels f1 and f2 of Figure 13). An analogous strain rate dependence of deformation was previously noted by Olive et al. (2016) in rifting simulations that treat the upper crust as a visco-plastic medium rather than an elasto-plastic medium. In the present case, however, the viscosity of the upper crust has a physical meaning ($\eta_D$, related to the damage growth rate), as opposed to an arbitrarily high value meant to simulate a stiff visco-plastic upper crust.

<figure><img src="image_17.png" /><figcaption>**Figure 15.** Numerical simulation of about 3.5 km of extension of a 10-km thick brittle plate, at varying rates (columns). Snapshots of the second invariant of strain rate at different amounts of finite extension (shown on the left), for a brittle plate governed by a–c the SCAM flow law, and d–f the EP-tecto rheology. The model resolution has been reduced compared to Figure 13 to reach greater amounts of finite extension in reasonable time.</figcaption></figure>
Snapshots of the second invariant of strain rate ($\dot{\varepsilon}_{II}$) fields at different times are displayed in Figure 15. The first three rows correspond to successive times throughout the SCAM simulations after ~251, ~1,079, and ~3,457 m of horizontal stretching, respectively. The fourth row corresponds to snapshots of elasto-plastic simulations performed under the same extension rates, after ~3,457 m of extension. The elasto-plastic simulations remain unsurprinsingly insensitive to a change of strain rate (e.g., Olive et al., 2016), as they all develop the same central graben structure (Figures 15d–15f).

The SCAM simulations, on the other hand, generate features that vary with extension rate. The first one is the location of the cluster of faults in the hanging-wall (Figure 15a1–15c1). The faster the extension rate, the closer this cluster of secondary faults lies to the master fault, with distances ranging from 15 to 7 km. The second feature is the geometry of the rift after 3,457 m of extension (Panels a3–c3). The rift structure shows a greater degree of asymmetry with increasing extension rate. This is particularly well expressed in topography, which corresponds to a symmetric graben in Panel a3 but is closer to a half-graben in Panels b3 and c3, as more extension is accommodated of the central right-dipping fault. We acknowledge that some of this variability could be attributed to stochasticity in marker positions, coupled with strong non-linearities of the SCAM flow law. We however ran the SCAM simulations a second time using different randomly assigned marker positions, and otherwise identical parameters. This second set yielded the same pattern of strain rate dependence as illustrated in Panels a3–c3.

Extension rate is known to influence rifting styles, primarily by modulating the thermal structure (Buck, 1991; Lavier & Buck, 2002). Our results suggest that it may play an additional, more subtle role by modulating the very processes of fault development. Specifically, the SCAM model introduces a strain rate dependence of brittle deformation by transiently activating a moderate to low viscosity in portions of the upper crust where damage is actively growing (e.g., Panels f1 and f2 of Figure 13). An analogous strain rate dependence of deformation was previously noted by Olive et al. (2016) in rifting simulations that treat the upper crust as a visco-plastic medium rather than an elasto-plastic medium. In the present case, however, the viscosity of the upper crust has a physical meaning ($\eta_D$, related to the damage growth rate), as opposed to an arbitrarily high value meant to simulate a stiff visco-plastic upper crust.

<figure><img src="image_17.png" /><figcaption>**Figure 15.** Numerical simulation of about 3.5 km of extension of a 10-km thick brittle plate, at varying rates (columns). Snapshots of the second invariant of strain rate at different amounts of finite extension (shown on the left), for a brittle plate governed by a–c the SCAM flow law, and d–f the EP-tecto rheology. The model resolution has been reduced compared to Figure 13 to reach greater amounts of finite extension in reasonable time.</figcaption></figure>
## 6.5. Upscaling the SCAM Parameters

The results presented thus far show that a flow law indexed on the activity of sub-millimetric rock defects, and calibrated on decimeter-sized samples can produce reasonable deformation patterns when applied to tectonic problems at scales >10 km. This fact does not negate the implication of larger defects (>10 cm) in the nucleation of crustal faults, but highlights that they need not be explicitly described in our constitutive relation. These intermediate-scale defects can be thought of as emerging through the formation of distributed damage bands, albeit only at scales greater than the numerical grid size, before dominant crustal-scale bands fully develop (e.g., Figure 13d1–13f1). In other words, the coalescence of sub-millimetric shear cracks produce diffuse proto-faults, and stress concentrations at their tips generate further localization at larger and larger scales.

Nonetheless, it is legitimate to wonder how differently our model would behave if it was based on larger (e.g., metric) shear defects. To address this question, we assess the effect of initial defect size *a* on the two extreme measures of rock strength in the SCAM framework: the minimum strength that characterizes extremely slow deformation, and the maximum strength at very high strain rates. The absolute lower bound on SCAM strength is given by the maximum differential stress satisfying $K_I = 0$ (Figure 6, end of Section 5.1.2). It corresponds to the differential stress that must be reached for macroscopic failure to be achievable, albeit after an infinitely long time (plain blue line in Figure 11). The expression of $K_I$ (19) is a function of $\sigma_1$, $\sigma_3$, $\mu$, $a$, and $N_v$ (through $D_0$), but in the special case where $K_I = 0$, the explicit dependence on $a$ drops, such that the maximum differential stress at $K_I = 0$ only depends on $\sigma_3$, $D_0$ and $\mu$. Assuming $\mu$ is constant across scales, we end up with a lower bound on strength that only depends on $D_0 = f(a, N_v)$.

Interestingly, this property can be used to constrain a plausible range of $D_0$ values to be used in tectonic simulations. The minimum strength (continuous blue line in Figure 11) must exceed the stress required to slip on a favorably oriented pre-existing frictional surfaces (e.g., ~500 MPa for $\sigma_3 = 150$ MPa, the value used to construct Figure 6; dashed blue line in Figure 11). The minimum strength must also be lower than the strength measured in the laboratory at the slowest possible strain rate (somewhere within the red area in Figure 11). Taking 90% of the peak value from Wawersik and Brace (1971) (e.g., 1,150 MPa for $\sigma_3 = 150$ MPa) as a rather conservative estimate, considering that Brantut et al. (2012) observed brittle creep failure at ~77% of estimated peak strength under 50 MPa of confining pressure, brackets the minimum strength for $\sigma_3 = 150$ MPa between ~500 and ~1,150 MPa. This corresponds to a range of $D_0$ between 0.06 and 0.4.

On the other hand, the maximum strength of the material at very fast strain rates can be assessed by equating Equation 19 to $K_{IC}$, and rewriting it as a function of $\sigma_1$. This gives an equation whose maximum with respect to *D* yields the maximum differential stress a rock can withstand (Figure 6), as cracks grow critically at seismic wave speeds (Bhat et al., 2012). This maximum stress roughly scales as $K_{IC}/\sqrt{a}$. It follows that the difference between the maximum differential stress at $K_I = K_{IC}$ and at $K_I = 0$, a measure of the overall strain rate dependence of the SCAM model, approximately scales as $a^{-1/2}$. If fracture toughness is a scale-invariant, intrinsic property of the material, increasing the size of the initial shear cracks brings the maximum and minimum strengths closer and closer to each other, severely suppressing the strain-rate dependence of the material's intact strength. This result is a fundamental property of the wing-crack model of Ashby and Sammis (1990), and does not depend on the specifics of the SCAM model. Whether it applies to natural systems remains an open question, as it is well known that fracture toughness can vary significantly across scales.

## 7. Conclusions and Perspectives

In this paper, we introduced a Sub-Critically Altered Maxwell (SCAM) framework to describe brittle deformation in long-term tectonic models. It is a set of constitutive equations that capture experimentally described behaviors of rocks at upper crustal pressures and temperatures. It is based on the evolution of an internal damage state and on its interactions with the elastic properties of the material. The model also allows large deformations by branching to plastic behavior after microcracks interact and coalesce. The SCAM model has several notable properties that make it a promising alternative to standard elasto-plastic models, or a way to improve their parameterizations of brittle yielding. Elastic properties are permanently altered due to their indexation on damage. Damage growth is a time-dependent process that is activated at stresses far below failure strength, through frictional sliding on shear defects distributed throughout the rock. It promotes strain softening pre-peak and self-consistently generates the successive stages of brittle creep. Creep behavior results from the transition from negative to positive retro-action
## 6.5. Upscaling the SCAM Parameters

The results presented thus far show that a flow law indexed on the activity of sub-millimetric rock defects, and calibrated on decimeter-sized samples can produce reasonable deformation patterns when applied to tectonic problems at scales >10 km. This fact does not negate the implication of larger defects (>10 cm) in the nucleation of crustal faults, but highlights that they need not be explicitly described in our constitutive relation. These intermediate-scale defects can be thought of as emerging through the formation of distributed damage bands, albeit only at scales greater than the numerical grid size, before dominant crustal-scale bands fully develop (e.g., Figure 13d1–13f1). In other words, the coalescence of sub-millimetric shear cracks produce diffuse proto-faults, and stress concentrations at their tips generate further localization at larger and larger scales.

Nonetheless, it is legitimate to wonder how differently our model would behave if it was based on larger (e.g., metric) shear defects. To address this question, we assess the effect of initial defect size *a* on the two extreme measures of rock strength in the SCAM framework: the minimum strength that characterizes extremely slow deformation, and the maximum strength at very high strain rates. The absolute lower bound on SCAM strength is given by the maximum differential stress satisfying $K_I = 0$ (Figure 6, end of Section 5.1.2). It corresponds to the differential stress that must be reached for macroscopic failure to be achievable, albeit after an infinitely long time (plain blue line in Figure 11). The expression of $K_I$ (19) is a function of $\sigma_1$, $\sigma_3$, $\mu$, $a$, and $N_v$ (through $D_0$), but in the special case where $K_I = 0$, the explicit dependence on $a$ drops, such that the maximum differential stress at $K_I = 0$ only depends on $\sigma_3$, $D_0$ and $\mu$. Assuming $\mu$ is constant across scales, we end up with a lower bound on strength that only depends on $D_0 = f(a, N_v)$.

Interestingly, this property can be used to constrain a plausible range of $D_0$ values to be used in tectonic simulations. The minimum strength (continuous blue line in Figure 11) must exceed the stress required to slip on a favorably oriented pre-existing frictional surfaces (e.g., ~500 MPa for $\sigma_3 = 150$ MPa, the value used to construct Figure 6; dashed blue line in Figure 11). The minimum strength must also be lower than the strength measured in the laboratory at the slowest possible strain rate (somewhere within the red area in Figure 11). Taking 90% of the peak value from Wawersik and Brace (1971) (e.g., 1,150 MPa for $\sigma_3 = 150$ MPa) as a rather conservative estimate, considering that Brantut et al. (2012) observed brittle creep failure at ~77% of estimated peak strength under 50 MPa of confining pressure, brackets the minimum strength for $\sigma_3 = 150$ MPa between ~500 and ~1,150 MPa. This corresponds to a range of $D_0$ between 0.06 and 0.4.

On the other hand, the maximum strength of the material at very fast strain rates can be assessed by equating Equation 19 to $K_{IC}$, and rewriting it as a function of $\sigma_1$. This gives an equation whose maximum with respect to *D* yields the maximum differential stress a rock can withstand (Figure 6), as cracks grow critically at seismic wave speeds (Bhat et al., 2012). This maximum stress roughly scales as $K_{IC}/\sqrt{a}$. It follows that the difference between the maximum differential stress at $K_I = K_{IC}$ and at $K_I = 0$, a measure of the overall strain rate dependence of the SCAM model, approximately scales as $a^{-1/2}$. If fracture toughness is a scale-invariant, intrinsic property of the material, increasing the size of the initial shear cracks brings the maximum and minimum strengths closer and closer to each other, severely suppressing the strain-rate dependence of the material's intact strength. This result is a fundamental property of the wing-crack model of Ashby and Sammis (1990), and does not depend on the specifics of the SCAM model. Whether it applies to natural systems remains an open question, as it is well known that fracture toughness can vary significantly across scales.

## 7. Conclusions and Perspectives

In this paper, we introduced a Sub-Critically Altered Maxwell (SCAM) framework to describe brittle deformation in long-term tectonic models. It is a set of constitutive equations that capture experimentally described behaviors of rocks at upper crustal pressures and temperatures. It is based on the evolution of an internal damage state and on its interactions with the elastic properties of the material. The model also allows large deformations by branching to plastic behavior after microcracks interact and coalesce. The SCAM model has several notable properties that make it a promising alternative to standard elasto-plastic models, or a way to improve their parameterizations of brittle yielding. Elastic properties are permanently altered due to their indexation on damage. Damage growth is a time-dependent process that is activated at stresses far below failure strength, through frictional sliding on shear defects distributed throughout the rock. It promotes strain softening pre-peak and self-consistently generates the successive stages of brittle creep. Creep behavior results from the transition from negative to positive retro-action
between damage and damage growth, which represents increasing interactions between lengthening cracks. The SCAM model also produces two yield strength using a single friction coefficient and no cohesion: a peak failure strength with high effective friction and strain rate dependent cohesion, which ultimately transitions to a rate-independent residual strength that obeys Byerlee's law. Despite the high effective friction (~1), shear band orientations in 2-D plane strain simulations remain consistent with the true friction coefficient that describes the strength of shear defects (μ ~ 0.7). The SCAM model can be calibrated against experimental data using prior knowledge on well-constrained rock properties to model the deformation of a specific lithology. Here Bayesian inversions of experimental data on Westerly granite and Darley Dale sandstone led to a set of reasonable micromechanical parameters whose impact on the macroscopic behavior of the rock can be straightforwardly interpreted.

Preliminary results of rifting simulations in a 10-km thick brittle plate subjected to gravity show that the SCAM model generates a population of faults with power-law distributed slip rates, akin to the distribution of natural fracture sizes and fault offsets, and likely introduces a strain rate dependence of the geological structures that develop at large strains. These features make it a good candidate to further investigate the complexity of brittle behavior across scales.

In designing the SCAM framework, we have strived to capture key micromechanical processes while keeping the model as simple as possible. This came at the cost of a few simplifications that may be relaxed in future work. A first strong assumption is that weakening of the shear modulus occurs isotropically, even though the orientation of the shear defects is strongly anisotropic. This assumption has the advantage of reducing the complexity of the constitutive law and simplifying its interpretation, but lacks consistency when relating crack geometry to macroscopic behavior. Within this assumption, however, a more consistent approach to damage-induced elastic weakening could be to assume isotropic weakening in the plane containing the crack normals (the ($σ_1, σ_3$) plane), and no weakening in the perpendicular direction. Such anisotropy however implies that the 2-D plane strain condition can no longer be satisfied for arbitrary extremum stresses, that is, enforcing $ε_2 = 0$ could lead to $σ_2$ not being an intermediate stress between $σ_1$ and $σ_3$.

In addition to neglecting material anisotropy, the evaluation of $K_I$ (Section 2.3) assumes that all cracks remain oriented at 45° with respect to $σ_1$ regardless of the stress and material rotations that may occur during large deformations. This assumption can lead to significant errors on the stresses resolved on shear defects. We note that the strains associated with shear band localization are on the order of a few percent, and thus should not be accompanied by significant finite rotation. On the other hand, large strains could conceivably rotate shear defects out of an optimal orientation in a long-term tectonic simulation. From a theoretical and computational point of view, keeping track of multiple shear defect orientations in an evolving stress state is challenging, in part because evaluating stress intensity factors and interaction terms would become very complex. The very definition of damage should also be revised to account for curving wing-cracks that grow along a changing $σ_1$ direction, and/or different sets of tensile wings associated with different sets of shear defects.

Another possible improvement for long-term tectonic modeling would be to allow damage to partially or completely heal over long time scales. This would promote the progressive recovery of intact elastic properties within abandoned shear bands, and require a re-mobilization of micro-scale frictional processes in order to re-activate a de-activated fault. To this end, a sub-critical crack growth law similar to that of Darot and Gueguen (1986) seems adequate. Crack growth is formulated as a thermally activated process where the growth or healing of a crack is indexed on the sign of the energy balance associated to the incremental advance of the crack front. This formulation has a more robust thermodynamical foundation than Charles' (1958) law and also accounts for the effect of temperature. It however involves more parameters, which is why we restricted ourselves to a simpler version of Charles' law in this work.

Finally, and importantly, the present study neglected the intrinsically dilatant effect of mode-I crack growth. Incompressibility is a common simplification in long-term tectonic models (Gerya, 2010), which amounts to setting the Poisson's ratio to 0.5 and the dilatancy angle to zero when modeling brittle failure with Mohr-Coulomb plasticity. Recent studies have shown that accounting for elastic and plastic compressibility noticeably impacts the outcome of tectonic simulations. Specifically, Duretz et al. (2021) found that compressibility hinders strain localization (it produces broader shear zones), and can facilitate the convergence of numerical solvers. A mechanically consistent way of introducing inelastic dilatancy within the elastically incompressible ($ν = 0.5$) SCAM model would be to use the constitutive law developed by Bhat et al. (2012). Their approach consisted of assessing
between damage and damage growth, which represents increasing interactions between lengthening cracks. The SCAM model also produces two yield strength using a single friction coefficient and no cohesion: a peak failure strength with high effective friction and strain rate dependent cohesion, which ultimately transitions to a rate-independent residual strength that obeys Byerlee's law. Despite the high effective friction (~1), shear band orientations in 2-D plane strain simulations remain consistent with the true friction coefficient that describes the strength of shear defects (μ ~ 0.7). The SCAM model can be calibrated against experimental data using prior knowledge on well-constrained rock properties to model the deformation of a specific lithology. Here Bayesian inversions of experimental data on Westerly granite and Darley Dale sandstone led to a set of reasonable micromechanical parameters whose impact on the macroscopic behavior of the rock can be straightforwardly interpreted.

Preliminary results of rifting simulations in a 10-km thick brittle plate subjected to gravity show that the SCAM model generates a population of faults with power-law distributed slip rates, akin to the distribution of natural fracture sizes and fault offsets, and likely introduces a strain rate dependence of the geological structures that develop at large strains. These features make it a good candidate to further investigate the complexity of brittle behavior across scales.

In designing the SCAM framework, we have strived to capture key micromechanical processes while keeping the model as simple as possible. This came at the cost of a few simplifications that may be relaxed in future work. A first strong assumption is that weakening of the shear modulus occurs isotropically, even though the orientation of the shear defects is strongly anisotropic. This assumption has the advantage of reducing the complexity of the constitutive law and simplifying its interpretation, but lacks consistency when relating crack geometry to macroscopic behavior. Within this assumption, however, a more consistent approach to damage-induced elastic weakening could be to assume isotropic weakening in the plane containing the crack normals (the ($σ_1, σ_3$) plane), and no weakening in the perpendicular direction. Such anisotropy however implies that the 2-D plane strain condition can no longer be satisfied for arbitrary extremum stresses, that is, enforcing $ε_2 = 0$ could lead to $σ_2$ not being an intermediate stress between $σ_1$ and $σ_3$.

In addition to neglecting material anisotropy, the evaluation of $K_I$ (Section 2.3) assumes that all cracks remain oriented at 45° with respect to $σ_1$ regardless of the stress and material rotations that may occur during large deformations. This assumption can lead to significant errors on the stresses resolved on shear defects. We note that the strains associated with shear band localization are on the order of a few percent, and thus should not be accompanied by significant finite rotation. On the other hand, large strains could conceivably rotate shear defects out of an optimal orientation in a long-term tectonic simulation. From a theoretical and computational point of view, keeping track of multiple shear defect orientations in an evolving stress state is challenging, in part because evaluating stress intensity factors and interaction terms would become very complex. The very definition of damage should also be revised to account for curving wing-cracks that grow along a changing $σ_1$ direction, and/or different sets of tensile wings associated with different sets of shear defects.

Another possible improvement for long-term tectonic modeling would be to allow damage to partially or completely heal over long time scales. This would promote the progressive recovery of intact elastic properties within abandoned shear bands, and require a re-mobilization of micro-scale frictional processes in order to re-activate a de-activated fault. To this end, a sub-critical crack growth law similar to that of Darot and Guéguen (1986) seems adequate. Crack growth is formulated as a thermally activated process where the growth or healing of a crack is indexed on the sign of the energy balance associated to the incremental advance of the crack front. This formulation has a more robust thermodynamical foundation than Charles' (1958) law and also accounts for the effect of temperature. It however involves more parameters, which is why we restricted ourselves to a simpler version of Charles' law in this work.

Finally, and importantly, the present study neglected the intrinsically dilatant effect of mode-I crack growth. Incompressibility is a common simplification in long-term tectonic models (Gerya, 2010), which amounts to setting the Poisson's ratio to 0.5 and the dilatancy angle to zero when modeling brittle failure with Mohr-Coulomb plasticity. Recent studies have shown that accounting for elastic and plastic compressibility noticeably impacts the outcome of tectonic simulations. Specifically, Duretz et al. (2021) found that compressibility hinders strain localization (it produces broader shear zones), and can facilitate the convergence of numerical solvers. A mechanically consistent way of introducing inelastic dilatancy within the elastically incompressible ($ν = 0.5$) SCAM model would be to use the constitutive law developed by Bhat et al. (2012). Their approach consisted of assessing
the Gibbs free energy of a damaged solid (assuming the (Ashby & Sammis, 1990) micro-mechanical model), and deriving it with respect to stress to yield effective elastic compliances. This approach outlines two sources of volumetric strain during damage accumulation: one directly results from the elastic compressibility of the material, which is altered by damage, the other represents a coupling between shear and volumetric deformation (a direct consequence of mode-I cracks being wedged open by shear defects). An elastically incompressible, yet inelastically dilatant SCAM model could focus on capturing this shear-volumetric coupling through an effective dilatancy angle that can be directly related to damage and stress Bhat et al. (2012). Dilatancy should cease as soon as cracks coalesce ($D \sim 1$), and give way to standard incompressible Mohr-Coulomb plasticity.

In its current state, the SCAM framework already unveils novel prospects for tectonic modelers. Keeping track of a damage field enables new connections between crustal scale simulations and key observables. Damage fields could straightforwardly be converted to seismic velocity maps for comparison with crustal tomography data. The distribution of damage along the model's free surface could also be used as a proxy for the erodibility of rocks exposed to weathering, and help understand variability in erosion rates across tectonically active landscapes (Gallen et al., 2015; Molnar et al., 2007). Finally, damage accumulation is intrinsically linked to an increase in porosity within the rock, and fracture connectivity is a primary control on fluid pathways. The evolving damage field could therefore be used as a proxy for the permeability of the brittle crust (e.g., Perol and Bhat (2016)) enabling coupled models of progressive brittle failure and fluid flow. This would enable the consideration of both poro-elastic feedbacks between fluid pressure and deformation, as well as alterations of mechanical properties through fluid-rock reactions.

# Appendix A: Bayesian Inversion Procedure

## A1. Mathematical Description

As illustrated in Section 3, our “forward” problem consists of predicting a vector of data points **d**, for example, a time series of stress or brittle creep strain rates, given a set of rock parameters (Table 1) stored in a vector **m**. This problem synthetically writes **d** = **g**(**m**). By contrast, the inverse problem consists of finding the distribution of model parameters that best fits known experimental data **d**<sub>obs</sub> with associated uncertainty. We adopt a least-squares approach in which all probability distributions quantifying uncertainties are Gaussian. Uncorrelated data uncertainty can thus be represented as a Gaussian distribution centered on **d**<sub>obs</sub> with a diagonal covariance matrix **C**<sub>d</sub>. Similarly, a-priori knowledge on the distribution of model parameters can be assumed to follow a Gaussian distribution centered on **m**<sub>prior</sub>, with diagonal covariance matrix **C**<sub>m</sub>.

We seek the best-fitting model ($\tilde{\mathbf{m}}$) which minimizes the log-likelihood function: a measure of distance between $\mathbf{d}_{\text{obs}}$ and $\mathbf{d}$, weighted by data uncertainty and prior knowledge:

$$S(\mathbf{m}) = (\mathbf{g}(\mathbf{m}) - \mathbf{d}_{\text{obs}})^T \mathbf{C}_{\text{d}}^{-1} (\mathbf{g}(\mathbf{m}) - \mathbf{d}_{\text{obs}}) + (\mathbf{m} - \mathbf{m}_{\text{prior}})^T \mathbf{C}_{\text{m}}^{-1} (\mathbf{m} - \mathbf{m}_{\text{prior}}). \quad (A1)$$

The log-likelihood is then minimized using a Gauss-Newton iterative algorithm (Tarantola, 2005). The iterative scheme, initialized at $\mathbf{m}_0 = \mathbf{m}_{\text{prior}}$, reads:

$$\mathbf{m}_{n+1} = \mathbf{m}_n - \kappa_n (\mathbf{G}_n^T \mathbf{C}_D^{-1} \mathbf{G}_n + \mathbf{C}_m^{-1})^{-1} (\mathbf{G}_n^T \mathbf{C}_D^{-1} (\mathbf{g}(\mathbf{m}_n) - \mathbf{d}_{\text{obs}}) + \mathbf{C}_m^{-1} (\mathbf{m}_n - \mathbf{m}_{\text{prior}})), \quad (A2)$$

where $n$ and $n + 1$ refer to the current and next iteration, $\kappa_n \le 1$ is the step multiplier, a hyperparameter that can be tuned to help convergence in case of strong non-linearities in the log-likelihood function. $\mathbf{G}$ is the Jacobian matrix:

$$G_{n}^{ij} = \left( \frac{\partial g_{i}}{\partial m_{j}} \right)_{\mathbf{m}_{n}}. \quad (A3)$$

The posterior model covariance matrix, a measure of the uncertainty on the inverted parameters, is then computed as
the Gibbs free energy of a damaged solid (assuming the (Ashby & Sammis, 1990) micro-mechanical model), and deriving it with respect to stress to yield effective elastic compliances. This approach outlines two sources of volumetric strain during damage accumulation: one directly results from the elastic compressibility of the material, which is altered by damage, the other represents a coupling between shear and volumetric deformation (a direct consequence of mode-I cracks being wedged open by shear defects). An elastically incompressible, yet inelastically dilatant SCAM model could focus on capturing this shear-volumetric coupling through an effective dilatancy angle that can be directly related to damage and stress Bhat et al. (2012). Dilatancy should cease as soon as cracks coalesce ($D \sim 1$), and give way to standard incompressible Mohr-Coulomb plasticity.

In its current state, the SCAM framework already unveils novel prospects for tectonic modelers. Keeping track of a damage field enables new connections between crustal scale simulations and key observables. Damage fields could straightforwardly be converted to seismic velocity maps for comparison with crustal tomography data. The distribution of damage along the model's free surface could also be used as a proxy for the erodibility of rocks exposed to weathering, and help understand variability in erosion rates across tectonically active landscapes (Gallen et al., 2015; Molnar et al., 2007). Finally, damage accumulation is intrinsically linked to an increase in porosity within the rock, and fracture connectivity is a primary control on fluid pathways. The evolving damage field could therefore be used as a proxy for the permeability of the brittle crust (e.g., Perol and Bhat (2016)) enabling coupled models of progressive brittle failure and fluid flow. This would enable the consideration of both poro-elastic feedbacks between fluid pressure and deformation, as well as alterations of mechanical properties through fluid-rock reactions.

# Appendix A: Bayesian Inversion Procedure

## A1. Mathematical Description

As illustrated in Section 3, our “forward” problem consists of predicting a vector of data points **d**, for example, a time series of stress or brittle creep strain rates, given a set of rock parameters (Table 1) stored in a vector **m**. This problem synthetically writes **d** = **g**(**m**). By contrast, the inverse problem consists of finding the distribution of model parameters that best fits known experimental data **d**<sub>obs</sub> with associated uncertainty. We adopt a least-squares approach in which all probability distributions quantifying uncertainties are Gaussian. Uncorrelated data uncertainty can thus be represented as a Gaussian distribution centered on **d**<sub>obs</sub> with a diagonal covariance matrix **C**<sub>d</sub>. Similarly, a-priori knowledge on the distribution of model parameters can be assumed to follow a Gaussian distribution centered on **m**<sub>prior</sub>, with diagonal covariance matrix **C**<sub>m</sub>.

We seek the best-fitting model ($\tilde{\mathbf{m}}$) which minimizes the log-likelihood function: a measure of distance between $\mathbf{d}_{\text{obs}}$ and $\mathbf{d}$, weighted by data uncertainty and prior knowledge:

$$S(\mathbf{m}) = (\mathbf{g}(\mathbf{m}) - \mathbf{d}_{\text{obs}})^T \mathbf{C}_{\text{d}}^{-1} (\mathbf{g}(\mathbf{m}) - \mathbf{d}_{\text{obs}}) + (\mathbf{m} - \mathbf{m}_{\text{prior}})^T \mathbf{C}_{\text{m}}^{-1} (\mathbf{m} - \mathbf{m}_{\text{prior}}). \quad (\text{A1})$$

The log-likelihood is then minimized using a Gauss-Newton iterative algorithm (Tarantola, 2005). The iterative scheme, initialized at $\mathbf{m}_0 = \mathbf{m}_{\text{prior}}$, reads:

$$\mathbf{m}_{n+1} = \mathbf{m}_n - \kappa_n (\mathbf{G}_n^T \mathbf{C}_D^{-1} \mathbf{G}_n + \mathbf{C}_m^{-1})^{-1} (\mathbf{G}_n^T \mathbf{C}_D^{-1} (\mathbf{g}(\mathbf{m}_n) - \mathbf{d}_{\text{obs}}) + \mathbf{C}_m^{-1} (\mathbf{m}_n - \mathbf{m}_{\text{prior}})), \quad (\text{A2})$$

where $n$ and $n + 1$ refer to the current and next iteration, $\kappa_n \le 1$ is the step multiplier, a hyperparameter that can be tuned to help convergence in case of strong non-linearities in the log-likelihood function. **G** is the Jacobian matrix:

$$G_{n}^{ij} = \left( \frac{\partial g_{i}}{\partial m_{j}} \right)_{\mathbf{m}_{n}}. \quad (\text{A3})$$

The posterior model covariance matrix, a measure of the uncertainty on the inverted parameters, is then computed as
$$ \tilde{C}_M \approx (G' C_D^{-1} G + C_M^{-1})^{-1} = C_M - C_M G' (G C_M G' + C_D)^{-1} G C_M. \quad (A4) $$

## A2. Implementation

For each inversion, we assemble a data vector $d_{\text{obs}}$ by concatenating data from several constant strain rate experiments conducted under different confining pressures, and several brittle creep experiments conducted under different imposed axial stresses. The constant strain rate data consists of time series of axial stress complemented by the peak axial stress and the corresponding axial strain that was reached in each experiment. By assigning a lower uncertainty (in $C_d$) on the peak axial stress and strain relative to the uncertainty on the stress time series, we can assign more weight to this constraint. Doing so helps favor models that accurately predict the position of the peak stress. The brittle creep data for Darley Dale sandstone consist of concatenated time series of axial strain rate from experiments conducted under different axial stresses. Because such data was not available for Westerly granite, we instead concatenated measurements of the representative secondary creep strain rate from experiments conducted under different axial stresses.

At every step of the inversion algorithm, a vector of “simulated data” $d_n = g(m_n)$ is built by simulating each individual experiment, and ordering the outputs to match the order of $d_{\text{obs}}$. In order to simulate brittle creep, the axial stress is first raised by imposing a constant strain rate matching that of the experiment (e.g., loading up to the star in Figure 4a as damage accumulates in Figure 4b). The damage state achieved by the sample is then used as initial condition for the brittle creep simulation. The axial stress is kept constant as we simulate the evolution of axial strain rates (See Section 3).

To ensure that the inversion does not assign negative values to parameters which are inherently positive, we construct a model vector $m$ that contains the logarithm of each SCAM parameter, namely: $G_0$, $\gamma$, $\mu$, $a$, $D_0$, $D_i$, $n$, $K_{IC}$, and $\dot{l}_0$. The associated uncertainties (prior or posterior) are thus log-normal distributions of the parameters, or Gaussian distributions of the logarithm of the parameters. If a positive parameter has a standard deviation $\bar{s}$ and a median $\bar{m}$, the distribution of its logarithm can be adequately represented by a Gaussian distribution centered on $\log(\bar{m})$ with variance $\bar{\sigma}^2 = \log(1 + \bar{s}^2/\bar{m}^2)$. We make use of these formula to assign uncertainties on the experimental data and on our prior knowledge of the model parameters.

## Appendix B: Derivation of the Plastic Viscosity

Two situations have to be considered to construct the plastic viscosity term $\eta_p$ used in SCAM numerical simulations. The first situation is stresses that lie above the plastic yield envelope at the onset of plastic behavior. This can happen because damage build-up may cause stresses to increase past the plastic yield stress ($\sigma_y$). In this case, we require the plastic viscosity to relax excess stresses above $\sigma_y$ in one iteration. Equation 38 recast in term of $\dot{s}_{ij}$ thus becomes

$$ \dot{s}_{ij} = \frac{s_{ij}^{y} - s_{ij}}{\Delta t} = 2G_{0}f(D)\left(\dot{e}_{ij} - \frac{s_{ij}}{2\eta_{p}}\right), \quad (B1) $$

where $s_{ij}^y$ corresponds to any deviatoric stress tensor whose second invariant satisfies $s_{II}^y = \sigma_y$. We approximate Equation B1 by using $\sigma_y$, $s_{II}$ and the second invariant of the strain rate tensors $\dot{e}_{II} = \sqrt{J_2(\dot{e})}$:

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

This yields the following closed form equation for $\eta_p$:

$$ \eta_p = \frac{s_{II}}{2\left(\dot{e}_{II} - \frac{\sigma_y - s_{II}}{2G_0 f(D) \Delta t}\right)}. \quad (B3) $$
$$ \tilde{C}_M \approx (G' C_D^{-1} G + C_M^{-1})^{-1} = C_M - C_M G' (G C_M G' + C_D)^{-1} G C_M. \quad (A4) $$

## A2. Implementation

For each inversion, we assemble a data vector $d_{\text{obs}}$ by concatenating data from several constant strain rate experiments conducted under different confining pressures, and several brittle creep experiments conducted under different imposed axial stresses. The constant strain rate data consists of time series of axial stress complemented by the peak axial stress and the corresponding axial strain that was reached in each experiment. By assigning a lower uncertainty (in $C_d$) on the peak axial stress and strain relative to the uncertainty on the stress time series, we can assign more weight to this constraint. Doing so helps favor models that accurately predict the position of the peak stress. The brittle creep data for Darley Dale sandstone consist of concatenated time series of axial strain rate from experiments conducted under different axial stresses. Because such data was not available for Westerly granite, we instead concatenated measurements of the representative secondary creep strain rate from experiments conducted under different axial stresses.

At every step of the inversion algorithm, a vector of “simulated data” $d_n = g(m_n)$ is built by simulating each individual experiment, and ordering the outputs to match the order of $d_{\text{obs}}$. In order to simulate brittle creep, the axial stress is first raised by imposing a constant strain rate matching that of the experiment (e.g., loading up to the star in Figure 4a as damage accumulates in Figure 4b). The damage state achieved by the sample is then used as initial condition for the brittle creep simulation. The axial stress is kept constant as we simulate the evolution of axial strain rates (See Section 3).

To ensure that the inversion does not assign negative values to parameters which are inherently positive, we construct a model vector $m$ that contains the logarithm of each SCAM parameter, namely: $G_0$, $\gamma$, $\mu$, $a$, $D_0$, $D_i$, $n$, $K_{IC}$, and $\dot{l}_0$. The associated uncertainties (prior or posterior) are thus log-normal distributions of the parameters, or Gaussian distributions of the logarithm of the parameters. If a positive parameter has a standard deviation $\bar{s}$ and a median $\bar{m}$, the distribution of its logarithm can be adequately represented by a Gaussian distribution centered on $\log(\bar{m})$ with variance $\bar{\sigma}^2 = \log(1 + \bar{s}^2/\bar{m}^2)$. We make use of these formula to assign uncertainties on the experimental data and on our prior knowledge of the model parameters.

## Appendix B: Derivation of the Plastic Viscosity

Two situations have to be considered to construct the plastic viscosity term $\eta_p$ used in SCAM numerical simulations. The first situation is stresses that lie above the plastic yield envelope at the onset of plastic behavior. This can happen because damage build-up may cause stresses to increase past the plastic yield stress ($\sigma_y$). In this case, we require the plastic viscosity to relax excess stresses above $\sigma_y$ in one iteration. Equation 38 recast in term of $\dot{s}_{ij}$ thus becomes

$$ \dot{s}_{ij} = \frac{s_{ij}^{y} - s_{ij}}{\Delta t} = 2G_{0}f(D)\left(\dot{e}_{ij} - \frac{s_{ij}}{2\eta_{p}}\right), \quad (B1) $$

where $s_{ij}^y$ corresponds to any deviatoric stress tensor whose second invariant satisfies $s_{II}^y = \sigma_y$. We approximate Equation B1 by using $\sigma_y$, $s_{II}$ and the second invariant of the strain rate tensors $\dot{e}_{II} = \sqrt{J_2(\dot{e})}$:

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

This yields the following closed form equation for $\eta_p$:

$$ \eta_p = \frac{s_{II}}{2\left(\dot{e}_{II} - \frac{\sigma_y - s_{II}}{2G_0 f(D) \Delta t}\right)}. \quad (B3) $$
The other possible situation involves stresses capped by the yield envelope (i.e., $s_{II} < \sigma_{y}$). In this case, we additionally require that $\eta_{p}$ be infinite if the elastic stress rate $2G_{0}f(D)\dot{e}_{II}$ is lower than the value required to reach the yield stress in one iteration ($\frac{\sigma_{y}-s_{II}}{\Delta t}$). This amounts to rewriting $\eta_{p}$:

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

to ensure that the plastic viscosity is finite only when $s_{II} \ge \sigma_{y}$, or if an elastic stress increment suddenly brings $s_{II}$ above $\sigma_{y}$.

## Notation

### Mohr-Coulomb plasticity

* **μ** friction coefficient

* **φ** (= arctan μ) friction angle on shear defects ($\phi_{m}$ at macroscopic scale)

* **C**<sub>m</sub> (macroscopic) cohesion

* **σ**<sub>y</sub> yield stress

* **σ**<sub>y</sub><sup>(max)</sup> intact plastic yield stress (determined by μ<sub>max</sub> and C<sub>max</sub>)

* **μ**<sub>max</sub> initial friction coefficient (in strain weakened Mohr-Coulomb plasticity)

* **C**<sub>max</sub> initial cohesion (in strain weakened Mohr-Coulomb plasticity)

* **σ**<sub>y</sub><sup>(min)</sup> fully weakened plastic yield stress (determined by μ<sub>min</sub> and C<sub>min</sub>)

* **μ**<sub>min</sub> fully weakened friction coefficient

* **C**<sub>min</sub> fully weakened cohesion

* **Δe**<sub>II</sub><sup>p</sup> accumulated plastic strain needed to fully weaken the frictional properties

### Damage mechanics

* **D** ∈ [D<sub>0</sub>, 1] damage internal state variable

* **D**<sub>0</sub> ∈ [0, 1] damage value corresponding to no tensile defect in the rock

* **D**<sub>i</sub> initial damage

* **D**<sub>c</sub> critical damage at the transition between the isolated crack regime and the interacting crack regime

* **γ** = f(D = 1) ratio of residual over reference shear modulus

* **N**<sub>v</sub> number of shear defects per unit volume

* **V**<sub>c</sub> characteristic volume per crack (1/N<sub>v</sub>)

* **A**<sub>c</sub> characteristic area per crack

* **A**<sub>b</sub> average area that separates neighboring cracks (bulk area in the (σ<sub>1</sub>, σ<sub>2</sub>) plane)

* **ψ** shear defect angle with respect to σ<sub>1</sub>

* **α** cos ψ

* **a** shear defect radius

* **l** tensile “wing” crack length

* **K**<sub>I</sub> mode I stress intensity factor
The other possible situation involves stresses capped by the yield envelope (i.e., $s_{II} < \sigma_{y}$). In this case, we additionally require that $\eta_{p}$ be infinite if the elastic stress rate $2G_{0}f(D)\dot{e}_{II}$ is lower than the value required to reach the yield stress in one iteration ($\frac{\sigma_{y}-s_{II}}{\Delta t}$). This amounts to rewriting $\eta_{p}$:

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

to ensure that the plastic viscosity is finite only when $s_{II} \ge \sigma_{y}$, or if an elastic stress increment suddenly brings $s_{II}$ above $\sigma_{y}$.

## Notation

### Mohr-Coulomb plasticity

* **μ** friction coefficient

* **φ** (= arctan μ) friction angle on shear defects ($\phi_{m}$ at macroscopic scale)

* **C**<sub>m</sub> (macroscopic) cohesion

* **σ**<sub>y</sub> yield stress

* **σ**<sub>y</sub><sup>(max)</sup> intact plastic yield stress (determined by μ<sub>max</sub> and C<sub>max</sub>)

* **μ**<sub>max</sub> initial friction coefficient (in strain weakened Mohr-Coulomb plasticity)

* **C**<sub>max</sub> initial cohesion (in strain weakened Mohr-Coulomb plasticity)

* **σ**<sub>y</sub><sup>(min)</sup> fully weakened plastic yield stress (determined by μ<sub>min</sub> and C<sub>min</sub>)

* **μ**<sub>min</sub> fully weakened friction coefficient

* **C**<sub>min</sub> fully weakened cohesion

* **Δe**<sub>II</sub><sup>p</sup> accumulated plastic strain needed to fully weaken the frictional properties

### Damage mechanics

* **D** ∈ [D<sub>0</sub>, 1] damage internal state variable

* **D**<sub>0</sub> ∈ [0, 1] damage value corresponding to no tensile defect in the rock

* **D**<sub>i</sub> initial damage

* **D**<sub>c</sub> critical damage at the transition between the isolated crack regime and the interacting crack regime

* **γ** = f(D = 1) ratio of residual over reference shear modulus

* **N**<sub>v</sub> number of shear defects per unit volume

* **V**<sub>c</sub> characteristic volume per crack (1/N<sub>v</sub>)

* **A**<sub>c</sub> characteristic area per crack

* **A**<sub>b</sub> average area that separates neighboring cracks (bulk area in the (σ<sub>1</sub>, σ<sub>2</sub>) plane)

* **ψ** shear defect angle with respect to σ<sub>1</sub>

* **α** cos ψ

* **a** shear defect radius

* **l** tensile “wing” crack length

* **K**<sub>I</sub> mode I stress intensity factor
<table><tr><td>K<sub>I</sub><sup>(w)</sup></td><td>mode I stress intensity factor due to the wedging force F<sub>w</sub></td></tr><tr><td>K<sub>I</sub><sup>(σ<sub>3</sub>)</sup></td><td>mode I stress intensity factor due to σ<sub>3</sub></td></tr><tr><td>K<sub>I</sub><sup>(i)</sup></td><td>mode I stress intensity factor due to interactions between cracks</td></tr><tr><td>σ<sub>3</sub><sup>i</sup></td><td>internal stress acting in the direction of σ<sub>3</sub> resulting from cracks interaction</td></tr><tr><td>K<sub>IC</sub></td><td>mode I fracture toughness</td></tr><tr><td>l̇<sub>0</sub></td><td>characteristic crack growth rate</td></tr><tr><td>n</td><td>Charles law exponent (corrosion index)</td></tr><tr><td>β</td><td>geometric regularization factor</td></tr><tr><td>A<sub>1</sub>, A<sub>3</sub></td><td>constants that depend on friction and the orientation of shear defects</td></tr></table>

## Stresses and strains

<table><tr><td>ε<sub>ij</sub></td><td>strain tensor</td></tr><tr><td>σ<sub>1</sub></td><td>most compressive principal stress</td></tr><tr><td>σ<sub>3</sub></td><td>least compressive principal stress</td></tr><tr><td>Δσ</td><td>= σ<sub>3</sub> - σ<sub>1</sub>, differential stress</td></tr><tr><td>Δσ<sup>e</sup></td><td>differential stress at K<sub>I</sub> = 0 and D = D<sub>0</sub></td></tr><tr><td>Δσ<sub>bc</sub></td><td>Minimum brittle strength</td></tr><tr><td>P</td><td>= -σ<sub>kk</sub>/3 pressure</td></tr><tr><td>p<sub>c</sub></td><td>= -σ<sub>3</sub> confining pressure in experiments</td></tr><tr><td>e<sub>ij</sub></td><td>deviatoric strain tensor</td></tr><tr><td>s<sub>ij</sub></td><td>deviatoric stress tensor</td></tr><tr><td>e<sub>ax</sub></td><td>deviatoric axial strain</td></tr><tr><td>s<sub>ax</sub></td><td>deviatoric axial stress</td></tr><tr><td>J<sub>2</sub>(X)</td><td>= (dev(X<sub>ij</sub>)dev(X<sub>ij</sub>))/2 second invariant of the deviator of second order tensor X</td></tr><tr><td>s<sub>II</sub></td><td>= √J<sub>2</sub>(σ<sub>ij</sub>) scalar shear stress magnitude</td></tr><tr><td>s<sub>ij</sub><sup>y</sup></td><td>deviatoric stress tensor satisfying the Mohr-Coulomb yield criterion</td></tr><tr><td>e<sub>II</sub></td><td>= √J<sub>2</sub>(ε<sub>ij</sub>) scalar shear strain magnitude</td></tr><tr><td>ε<sub>ij</sub><sup>p</sup></td><td>plastic strain tensor</td></tr></table>

## Additional notations

<table><tr><td>G</td><td>effective shear modulus</td></tr><tr><td>G<sub>0</sub></td><td>= G(D = D<sub>0</sub>) reference shear modulus corresponding to the lowest damage state (no tensile defect)</td></tr><tr><td>f(D)</td><td>weakening function</td></tr><tr><td>ν</td><td>Poisson's ratio</td></tr><tr><td>η<sub>D</sub></td><td>damage viscosity</td></tr><tr><td>η<sub>p</sub></td><td>plastic viscosity (2-D SCAM simulations)</td></tr><tr><td>η<sub>eff</sub></td><td>effective viscosity (2-D SCAM simulations)</td></tr></table>
<table><tr><td>K<sub>I</sub><sup>(w)</sup></td><td>mode I stress intensity factor due to the wedging force F<sub>w</sub></td></tr><tr><td>K<sub>I</sub><sup>(σ<sub>3</sub>)</sup></td><td>mode I stress intensity factor due to σ<sub>3</sub></td></tr><tr><td>K<sub>I</sub><sup>(i)</sup></td><td>mode I stress intensity factor due to interactions between cracks</td></tr><tr><td>σ<sub>3</sub><sup>i</sup></td><td>internal stress acting in the direction of σ<sub>3</sub> resulting from cracks interaction</td></tr><tr><td>K<sub>IC</sub></td><td>mode I fracture toughness</td></tr><tr><td>l̇<sub>0</sub></td><td>characteristic crack growth rate</td></tr><tr><td>n</td><td>Charles law exponent (corrosion index)</td></tr><tr><td>β</td><td>geometric regularization factor</td></tr><tr><td>A<sub>1</sub>, A<sub>3</sub></td><td>constants that depend on friction and the orientation of shear defects</td></tr></table>

## Stresses and strains

<table><tr><td>ε<sub>ij</sub></td><td>strain tensor</td></tr><tr><td>σ<sub>1</sub></td><td>most compressive principal stress</td></tr><tr><td>σ<sub>3</sub></td><td>least compressive principal stress</td></tr><tr><td>Δσ</td><td>= σ<sub>3</sub> - σ<sub>1</sub>, differential stress</td></tr><tr><td>Δσ<sup>e</sup></td><td>differential stress at K<sub>I</sub> = 0 and D = D<sub>0</sub></td></tr><tr><td>Δσ<sub>bc</sub></td><td>Minimum brittle strength</td></tr><tr><td>P</td><td>= -σ<sub>kk</sub>/3 pressure</td></tr><tr><td>p<sub>c</sub></td><td>= -σ<sub>3</sub> confining pressure in experiments</td></tr><tr><td>e<sub>ij</sub></td><td>deviatoric strain tensor</td></tr><tr><td>s<sub>ij</sub></td><td>deviatoric stress tensor</td></tr><tr><td>e<sub>ax</sub></td><td>deviatoric axial strain</td></tr><tr><td>s<sub>ax</sub></td><td>deviatoric axial stress</td></tr><tr><td>J<sub>2</sub>(X)</td><td>= (dev(X<sub>ij</sub>)dev(X<sub>ij</sub>))/2 second invariant of the deviator of second order tensor X</td></tr><tr><td>s<sub>II</sub></td><td>= √J<sub>2</sub>(σ<sub>ij</sub>) scalar shear stress magnitude</td></tr><tr><td>s<sub>ij</sub><sup>y</sup></td><td>deviatoric stress tensor satisfying the Mohr-Coulomb yield criterion</td></tr><tr><td>e<sub>II</sub></td><td>= √J<sub>2</sub>(ε<sub>ij</sub>) scalar shear strain magnitude</td></tr><tr><td>ε<sub>ij</sub><sup>p</sup></td><td>plastic strain tensor</td></tr></table>

## Additional notations

<table><tr><td>G</td><td>effective shear modulus</td></tr><tr><td>G<sub>0</sub></td><td>= G(D = D<sub>0</sub>) reference shear modulus corresponding to the lowest damage state (no tensile defect)</td></tr><tr><td>f(D)</td><td>weakening function</td></tr><tr><td>ν</td><td>Poisson's ratio</td></tr><tr><td>η<sub>D</sub></td><td>damage viscosity</td></tr><tr><td>η<sub>p</sub></td><td>plastic viscosity (2-D SCAM simulations)</td></tr><tr><td>η<sub>eff</sub></td><td>effective viscosity (2-D SCAM simulations)</td></tr></table>
$\eta_{\min}$ minimum viscosity (2-D SCAM simulations)

$\eta_{\max}$ maximum viscosity (2-D SCAM simulations)

$v_i$ Components of the velocity field

$\rho$ density

$g_i$ Components of the gravity field

$\Delta t$ time step

## Data Availability Statement

The codes used to produce the results of Section 3 (0-D, written in Julia) and Sections 4 and 6 (2-D, written in MATLAB) are provided in a Zenodo repository (Petit & Olive, 2024).

## Acknowledgments

This work was partially supported by an Emergence(s) Ville de Paris grant to J.-A. O. H.S.B. acknowledges the European Research Council Grant PERSISMO (Grant 865411) for partial support of this work. We are grateful to Muriel Gerbault, Mike Heap, Cécile Prigent and Luc Lavier for their helpful comments on earlier versions of this work. We also thank our editor Boris Kaus as well as John Naliboff and an anonymous reviewer for their feedback.

## References

Aben, F. M., Brantut, N., Mitchell, T. M., & David, E. C. (2019). Rupture energetics in crustal rock from laboratory-scale seismic tomography. *Geophysical Research Letters*, 46(13), 7337–7344. https://doi.org/10.1029/2019GL083040

Arthur, J., Dunstan, T., Al-Ani, Q. A. J. L., & Assadi, A. (1977). Plastic deformation and failure in granular media. *Géotechnique*, 27(1), 53–74. https://doi.org/10.1680/geot.1977.27.1.53

Ashby, M. F., & Hallam, S. D. (1986). The failure of brittle solids containing small cracks under compressive stress states. *Acta Metallurgica*, 34(3), 497–510. https://doi.org/10.1016/0001-6160(86)90086-6

Ashby, M. F., & Sammis, C. G. (1990). The damage mechanics of brittle solids in compression. *Pure and Applied Geophysics* PAGEOPH, 133(3), 33–521. https://doi.org/10.1007/bf00878002

Atkinson, B. K. (1979). A fracture mechanics study of subcritical tensile cracking of quartz in wet environments. *Pure and Applied Geophysics* PAGEOPH, 117(5), 1011–1024. https://doi.org/10.1007/bf00876082

Atkinson, B. K. (1984). Subcritical crack growth in geological materials. *Journal of Geophysical Research*, 89(B6), 4077–4114. https://doi.org/10.1029/JB089iB06p04077

Baud, P., & Meredith, P. (1997). Damage accumulation during triaxial creep of Darley Dale sandstone from pore volumometry and acoustic emission. *International Journal of Rock Mechanics and Mining Sciences*, 34(3–4), 24–el. https://doi.org/10.1016/s1365-1609(97)00060-9

Baud, P., Schubnel, A., & Wong, T.-f. (2000). Dilatancy, compaction, and failure mode in Solnhofen limestone. *Journal of Geophysical Research*, 105(B8), 19289–19303. https://doi.org/10.1029/2000JB900133

Bhat, H. S., Rosakis, A. J., & Sammis, C. G. (2012). A micromechanics based constitutive model for brittle failure at high strain rates. *Journal of Applied Mechanics*, 79(3), 031016. https://doi.org/10.1115/1.4005897

Bhat, H. S., Sammis, C., & Rosakis, A. (2011). The micromechanics of westerley granite at large compressive loads. *Pure and Applied Geophysics*, 168(12), 2181–2198. https://doi.org/10.1007/s00024-011-0271-9

Brace, W. F., & Kohlstedt, D. L. (1980). Limits on lithospheric stress imposed by laboratory experiments. *Journal of Geophysical Research*, 85(B11), 6248–6252. https://doi.org/10.1029/JB085iB11p06248

Brace, W. F., Paulding, B. W., & Scholz, C. (1966). Dilatancy in the fracture of crystalline rocks. *Journal of Geophysical Research*, 71(16), 3939–3953. https://doi.org/10.1029/JZ071i016p03939

Brantut, N., Baud, P., Heap, M. J., & Meredith, P. G. (2012). Micromechanics of brittle creep in rocks: Micromechanics of brittle creep. *Journal of Geophysical Research*, 117(B8). https://doi.org/10.1029/2012JB009299

Brantut, N., Heap, M., Meredith, P., & Baud, P. (2013). Time-dependent cracking and brittle creep in crustal rocks: A review. *Journal of Structural Geology*, 52, 17–43. https://doi.org/10.1016/j.jsg.2013.03.007

Buck, W. R. (1991). Modes of continental lithospheric extension. *Journal of Geophysical Research*, 96(B12), 20161–20178. https://doi.org/10.1029/91JB01485

Budiansky, B., & O’connell, R. J. (1976). Elastic moduli of a cracked solid. *International Journal of Solids and Structures*, 12(2), 81–97. https://doi.org/10.1016/0020-7683(76)90044-5

Byerlee, J. D. (1967). Frictional characteristics of granite under high confining pressure. *Journal of Geophysical Research*, 72(14), 3639–3648. https://doi.org/10.1029/JZ072i014p03639

Byerlee, J. D. (1978). Friction of rocks. *Pure and Applied Geophysics* PAGEOPH, 116, 615–626.

Carter, N., Anderson, D. A., Hansen, F. D., & Kranz, R. L. (1981). Creep and creep rupture of granitic rocks. In *Washington DC American Geophysical Union Geophysical Monograph Series* (Vol. 24, pp. 61–82).

Charles, R. J. (1958). Static Fatigue of Glass. *I. Journal of Applied Physics*, 29(11), 1549–1553. https://doi.org/10.1063/1.1722991

Costin, L. (1983). A microcrack model for the deformation and failure of brittle rock. *Journal of Geophysical Research*, 88(B11), 9485–9492. https://doi.org/10.1029/JB088iB11p09485

Costin, L. (1985). Damage mechanics in the post-failure regime. *Mechanics of Materials*, 4(2), 149–160. https://doi.org/10.1016/0167-6636(85)90013-4

Darot, M., & Gueguen, Y. (1986). Slow crack growth in minerals and rocks: Theory and experiments. *Pure and Applied Geophysics* PAGEOPH, 124(4–5), 677–692. https://doi.org/10.1007/BF00879604

Deshpande, V., & Evans, A. (2008). Inelastic deformation and energy dissipation in ceramics: A mechanism-based constitutive model. *Journal of the Mechanics and Physics of Solids*, 56(10), 3077–3100. https://doi.org/10.1016/j.jmps.2008.05.002

Dey, T. N., & Wang, C.-Y. (1981). Some mechanisms of microcrack growth and interaction in compressive rock failure. *International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts*, 18(3), 199–209. https://doi.org/10.1016/0148-9062(81)90974-8
$\eta_{\min}$ minimum viscosity (2-D SCAM simulations)

$\eta_{\max}$ maximum viscosity (2-D SCAM simulations)

$v_i$ Components of the velocity field

$\rho$ density

$g_i$ Components of the gravity field

$\Delta t$ time step

## Data Availability Statement

The codes used to produce the results of Section 3 (0-D, written in Julia) and Sections 4 and 6 (2-D, written in MATLAB) are provided in a Zenodo repository (Petit & Olive, 2024).

## Acknowledgments

This work was partially supported by an Emergence(s) Ville de Paris grant to J.-A. O. H.S.B. acknowledges the European Research Council Grant PERSISMO (Grant 865411) for partial support of this work. We are grateful to Muriel Gerbault, Mike Heap, Cécile Prigent and Luc Lavier for their helpful comments on earlier versions of this work. We also thank our editor Boris Kaus as well as John Naliboff and an anonymous reviewer for their feedback.

## References

Aben, F. M., Brantut, N., Mitchell, T. M., & David, E. C. (2019). Rupture energetics in crustal rock from laboratory-scale seismic tomography. *Geophysical Research Letters*, 46(13), 7337–7344. https://doi.org/10.1029/2019GL083040

Arthur, J., Dunstan, T., Al-Ani, Q. A. J. L., & Assadi, A. (1977). Plastic deformation and failure in granular media. *Géotechnique*, 27(1), 53–74. https://doi.org/10.1680/geot.1977.27.1.53

Ashby, M. F., & Hallam, S. D. (1986). The failure of brittle solids containing small cracks under compressive stress states. *Acta Metallurgica*, 34(3), 497–510. https://doi.org/10.1016/0001-6160(86)90086-6

Ashby, M. F., & Sammis, C. G. (1990). The damage mechanics of brittle solids in compression. *Pure and Applied Geophysics* PAGEOPH, 133(3), 33–521. https://doi.org/10.1007/bf00878002

Atkinson, B. K. (1979). A fracture mechanics study of subcritical tensile cracking of quartz in wet environments. *Pure and Applied Geophysics* PAGEOPH, 117(5), 1011–1024. https://doi.org/10.1007/bf00876082

Atkinson, B. K. (1984). Subcritical crack growth in geological materials. *Journal of Geophysical Research*, 89(B6), 4077–4114. https://doi.org/10.1029/JB089iB06p04077

Baud, P., & Meredith, P. (1997). Damage accumulation during triaxial creep of Darley Dale sandstone from pore volumometry and acoustic emission. *International Journal of Rock Mechanics and Mining Sciences*, 34(3–4), 24–el. https://doi.org/10.1016/s1365-1609(97)00060-9

Baud, P., Schubnel, A., & Wong, T.-f. (2000). Dilatancy, compaction, and failure mode in Solnhofen limestone. *Journal of Geophysical Research*, 105(B8), 19289–19303. https://doi.org/10.1029/2000JB900133

Bhat, H. S., Rosakis, A. J., & Sammis, C. G. (2012). A micromechanics based constitutive model for brittle failure at high strain rates. *Journal of Applied Mechanics*, 79(3), 031016. https://doi.org/10.1115/1.4005897

Bhat, H. S., Sammis, C., & Rosakis, A. (2011). The micromechanics of westerley granite at large compressive loads. *Pure and Applied Geophysics*, 168(12), 2181–2198. https://doi.org/10.1007/s00024-011-0271-9

Brace, W. F., & Kohlstedt, D. L. (1980). Limits on lithospheric stress imposed by laboratory experiments. *Journal of Geophysical Research*, 85(B11), 6248–6252. https://doi.org/10.1029/JB085iB11p06248

Brace, W. F., Paulding, B. W., & Scholz, C. (1966). Dilatancy in the fracture of crystalline rocks. *Journal of Geophysical Research*, 71(16), 3939–3953. https://doi.org/10.1029/JZ071i016p03939

Brantut, N., Baud, P., Heap, M. J., & Meredith, P. G. (2012). Micromechanics of brittle creep in rocks: Micromechanics of brittle creep. *Journal of Geophysical Research*, 117(B8). https://doi.org/10.1029/2012JB009299

Brantut, N., Heap, M., Meredith, P., & Baud, P. (2013). Time-dependent cracking and brittle creep in crustal rocks: A review. *Journal of Structural Geology*, 52, 17–43. https://doi.org/10.1016/j.jsg.2013.03.007

Buck, W. R. (1991). Modes of continental lithospheric extension. *Journal of Geophysical Research*, 96(B12), 20161–20178. https://doi.org/10.1029/91JB01485

Budiansky, B., & O’connell, R. J. (1976). Elastic moduli of a cracked solid. *International Journal of Solids and Structures*, 12(2), 81–97. https://doi.org/10.1016/0020-7683(76)90044-5

Byerlee, J. D. (1967). Frictional characteristics of granite under high confining pressure. *Journal of Geophysical Research*, 72(14), 3639–3648. https://doi.org/10.1029/JZ072i014p03639

Byerlee, J. D. (1978). Friction of rocks. *Pure and Applied Geophysics* PAGEOPH, 116, 615–626.

Carter, N., Anderson, D. A., Hansen, F. D., & Kranz, R. L. (1981). Creep and creep rupture of granitic rocks. In *Washington DC American Geophysical Union Geophysical Monograph Series* (Vol. 24, pp. 61–82).

Charles, R. J. (1958). Static Fatigue of Glass. *I. Journal of Applied Physics*, 29(11), 1549–1553. https://doi.org/10.1063/1.1722991

Costin, L. (1983). A microcrack model for the deformation and failure of brittle rock. *Journal of Geophysical Research*, 88(B11), 9485–9492. https://doi.org/10.1029/JB088iB11p09485

Costin, L. (1985). Damage mechanics in the post-failure regime. *Mechanics of Materials*, 4(2), 149–160. https://doi.org/10.1016/0167-6636(85)90013-4

Darot, M., & Gueguen, Y. (1986). Slow crack growth in minerals and rocks: Theory and experiments. *Pure and Applied Geophysics* PAGEOPH, 124(4–5), 677–692. https://doi.org/10.1007/BF00879604

Deshpande, V., & Evans, A. (2008). Inelastic deformation and energy dissipation in ceramics: A mechanism-based constitutive model. *Journal of the Mechanics and Physics of Solids*, 56(10), 3077–3100. https://doi.org/10.1016/j.jmps.2008.05.002

Dey, T. N., & Wang, C.-Y. (1981). Some mechanisms of microcrack growth and interaction in compressive rock failure. *International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts*, 18(3), 199–209. https://doi.org/10.1016/0148-9062(81)90974-8
Duretz, T., Borst, R., & Yamato, P. (2021). Modeling lithospheric deformation using a compressible Visco-Elasto-Viscoplastic rheology and the effective viscosity approach. *Geochemistry, Geophysics, Geosystems*, 22(8). https://doi.org/10.1029/2021GC009675

Eppes, M.-C., & Keanini, R. (2017). Mechanical weathering and rock erosion by climate-dependent subcritical cracking. *Reviews of Geophysics*, 55(2), 470–508. https://doi.org/10.1002/2017RG000557

Gallen, S. F., Clark, M. K., & Godt, J. W. (2015). Coseismic landslides reveal near-surface rock strength in a high-relief, tectonically active setting. *Geology*, 43(1), 11–14. https://doi.org/10.1130/G36080.1

Gerya, T. (2010). *Introduction to numerical geodynamic modelling*. Cambridge University Press.

Gerya, T., & Yuen, D. A. (2003). Characteristics-based marker-in-cell method with conservative finite-differences schemes for modeling geological flows with strongly variable transport properties. *Physics of the Earth and Planetary Interiors*, 140(4), 293–318. https://doi.org/10.1016/j.pepi.2003.09.006

Gratier, J.-P., Dysthe, D. K., & Renard, F. (2013). The role of pressure solution creep in the ductility of the Earth's upper crust. *Advances in Geophysics*, 54, 47–179. https://doi.org/10.1016/B978-0-12-380940-7.00002-0

Hamiel, Y., Liu, Y., Lyakhovsky, V., Ben-Zion, Y., & Lockner, D. (2004). A viscoelastic damage model with applications to stable and unstable fracturing. *Geophysical Journal International*, 159(3), 1155–1165. https://doi.org/10.1111/j.1365-246X.2004.02452.x

Heap, M. J., Baud, P., Meredith, P. G., Bell, A. F., & Main, I. G. (2009). Time-dependent brittle creep in Darley Dale sandstone. *Journal of Geophysical Research*, 114(B7), B07203. https://doi.org/10.1029/2008JB006212

Kachanov, M. L. (1982a). A microcrack model of rock inelasticity part I: Frictional sliding on microcracks. *Mechanics of Materials*, 1(1), 19–27. https://doi.org/10.1016/0167-6636(82)90021-7

Kachanov, M. L. (1982b). Microcrack model of rock inelasticity part III: Time-dependent growth of microcracks. *Mechanics of Materials*, 1(2), 123–129. https://doi.org/10.1016/0167-6636(82)90040-0

Kachanov, M. L. (1982c). A microcrack model of rock inelasticity part II: Propagation of microcracks. *Mechanics of Materials*, 1(1), 29–41. https://doi.org/10.1016/0167-6636(82)90022-9

Kachanov, M. L. (1993). Elastic solids with many cracks and related problems. *Advances in Applied Mechanics*, 30, 259–445. https://doi.org/10.1016/S0065-2156(08)70176-5

Karabulut, H., & Bouchon, M. (2007). Spatial variability and non-linearity of strong ground motion near a fault. *Geophysical Journal International*, 170(1), 262–274. https://doi.org/10.1111/j.1365-246X.2007.03406.x

Karrech, A., Regenauer-Lieb, K., & Poulet, T. (2011). A damaged visco-plasticity model for pressure and temperature sensitive geomaterials. *International Journal of Engineering Science*, 49(10), 1141–1150. https://doi.org/10.1016/j.ijengsci.2011.05.005

Kaus, B. J. (2010). Factors that control the angle of shear bands in geodynamic numerical models of brittle deformation. *Tectonophysics*, 484(1–4), 36–47. https://doi.org/10.1016/j.tecto.2009.08.042

Kranz, R. L. (1979). Crack growth and development during creep of Barre granite. *International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts*, 16(1), 23–35. https://doi.org/10.1016/0148-9062(79)90772-1

Lavier, L. L., & Buck, W. R. (2002). Half graben versus large-offset low-angle normal fault: Importance of keeping cool during normal faulting. *Journal of Geophysical Research*, 107(B6), ETG–8. https://doi.org/10.1029/2001jb000513

Lavier, L. L., Buck, W. R., & Poliakov, A. N. B. (2000). Factors controlling normal fault offset in an ideal brittle layer. *Journal of Geophysical Research*, 105(B10), 23431–23442. https://doi.org/10.1029/2000JB00108

Lemiale, V., Mühlhaus, H.-B., Moresi, L., & Stafford, J. (2008). Shear banding analysis of plastic models formulated for incompressible viscous flows. *Physics of the Earth and Planetary Interiors*, 171(1–4), 177–186. https://doi.org/10.1016/j.pepi.2008.07.038

Le Pourhiet, L. (2013). Strain localization due to structural softening during pressure sensitive rate independent yielding. *Bulletin de la Societe Geologique de France*, 184(4–5), 357–371. https://doi.org/10.2113/gsgfbull.184.4-5.357

Lockner, D. A. (1998). A generalized law for brittle deformation of Westerly granite. *Journal of Geophysical Research*, 103(B3), 5107–5123. https://doi.org/10.1029/97JB03211

Lockner, D. A., Byerlee, J. D., Kuksenko, V., Ponomarev, A., & Sidorin, A. (1991). Quasi-static fault growth and shear fracture energy in granite. *Nature*, 350(6313), 39–42. https://doi.org/10.1038/350039a0

Lyakhovsky, V., Ben-Zion, Y., & Agnon, A. (1997). Distributed damage, faulting, and friction. *Journal of Geophysical Research*, 102(B12), 27635–27649. https://doi.org/10.1029/97JB01896

Manaker, D. M., Turcotte, D. L., & Kellogg, L. H. (2006). Flexure with damage. *Geophysical Journal International*, 166(3), 1368–1383. https://doi.org/10.1111/j.1365-246X.2006.03067.x

Marrett, R., Ortega, O. J., & Kelsey, C. M. (1999). Extent of power-law scaling for natural fractures in rock. *Geology*, 27(9), 799. https://doi.org/10.1130/0091-7613(1999)027(0799:EOPLSF)2.3.CO;2

McBeck, J., Kandula, N., Aiken, J. M., Cordonnier, B., & Renard, F. (2019). Isolating the factors that govern fracture development in rocks throughout dynamic in situ X-ray tomography experiments. *Geophysical Research Letters*, 46(20), 11127–11135. https://doi.org/10.1029/2019GL084613

Menéndez, B., Zhu, W., & Wong, T.-F. (1996). Micromechanics of brittle faulting and cataclastic flow in Berea sandstone. *Journal of Structural Geology*, 18(1), 1–16. https://doi.org/10.1016/0191-8141(95)00076-P

Meyer, S. E., Kaus, B. J. P., & Passchier, C. (2017). Development of branching brittle and ductile shear zones: A numerical study. *Geochemistry, Geophysics, Geosystems*, 18(6), 2054–2075. https://doi.org/10.1002/2016GC006793

Molnar, P., Anderson, R. S., & Anderson, S. P. (2007). Tectonics, fracturing of rock, and erosion. *Journal of Geophysical Research*, 112(F3). https://doi.org/10.1029/2005JF000433

Moresi, L., Dufour, F., & Mühlhaus, H.-B. (2003). A Lagrangian integration point finite element method for large deformation modeling of viscoelastic geomaterials. *Journal of Computational Physics*, 184(2), 476–497. https://doi.org/10.1016/S0021-9991(02)00031-1

Morley, C. K. (1996). Discussion of potential errors in fault heave methods for extension estimates in rifts, with particular reference to fractal fault populations and inherited fabrics. *Geological Society, London, Special Publications*, 99(1), 117–134. https://doi.org/10.1144/GSL.SP.1996.099.01.10

Naliboff, J. B., Glerum, A., Brune, S., Péron-Pinvidic, G., & Wrona, T. (2020). Development of 3-d rift heterogeneity through fault network evolution. *Geophysical Research Letters*, 47(13), e2019GL086611. https://doi.org/10.1029/2019GL086611

Nemat-Nasser, S., & Horii, H. (1982). Compression-induced nonplanar crack extension with application to splitting, exfoliation, and rockburst. *Journal of Geophysical Research*, 87(B8), 6805–6821. https://doi.org/10.1029/JB087iB08p06805

Olive, J.-A., Behn, M. D., Mittelstaedt, E., Ito, G., & Klein, B. Z. (2016). The role of elasticity in simulating long-term tectonic extension. *Geophysical Journal International*, 205(2), 728–743. https://doi.org/10.1093/gji/ggw044
Duretz, T., Borst, R., & Yamato, P. (2021). Modeling lithospheric deformation using a compressible Visco-Elasto-Viscoplastic rheology and the effective viscosity approach. *Geochemistry, Geophysics, Geosystems*, 22(8). https://doi.org/10.1029/2021GC009675

Eppes, M.-C., & Keanini, R. (2017). Mechanical weathering and rock erosion by climate-dependent subcritical cracking. *Reviews of Geophysics*, 55(2), 470–508. https://doi.org/10.1002/2017RG000557

Gallen, S. F., Clark, M. K., & Godt, J. W. (2015). Coseismic landslides reveal near-surface rock strength in a high-relief, tectonically active setting. *Geology*, 43(1), 11–14. https://doi.org/10.1130/G36080.1

Gerya, T. (2010). *Introduction to numerical geodynamic modelling*. Cambridge University Press.

Gerya, T., & Yuen, D. A. (2003). Characteristics-based marker-in-cell method with conservative finite-differences schemes for modeling geological flows with strongly variable transport properties. *Physics of the Earth and Planetary Interiors*, 140(4), 293–318. https://doi.org/10.1016/j.pepi.2003.09.006

Gratier, J.-P., Dysthe, D. K., & Renard, F. (2013). The role of pressure solution creep in the ductility of the Earth's upper crust. *Advances in Geophysics*, 54, 47–179. https://doi.org/10.1016/B978-0-12-380940-7.00002-0

Hamiel, Y., Liu, Y., Lyakhovsky, V., Ben-Zion, Y., & Lockner, D. (2004). A viscoelastic damage model with applications to stable and unstable fracturing. *Geophysical Journal International*, 159(3), 1155–1165. https://doi.org/10.1111/j.1365-246X.2004.02452.x

Heap, M. J., Baud, P., Meredith, P. G., Bell, A. F., & Main, I. G. (2009). Time-dependent brittle creep in Darley Dale sandstone. *Journal of Geophysical Research*, 114(B7), B07203. https://doi.org/10.1029/2008JB006212

Kachanov, M. L. (1982a). A microcrack model of rock inelasticity part I: Frictional sliding on microcracks. *Mechanics of Materials*, 1(1), 19–27. https://doi.org/10.1016/0167-6636(82)90021-7

Kachanov, M. L. (1982b). Microcrack model of rock inelasticity part III: Time-dependent growth of microcracks. *Mechanics of Materials*, 1(2), 123–129. https://doi.org/10.1016/0167-6636(82)90040-0

Kachanov, M. L. (1982c). A microcrack model of rock inelasticity part II: Propagation of microcracks. *Mechanics of Materials*, 1(1), 29–41. https://doi.org/10.1016/0167-6636(82)90022-9

Kachanov, M. L. (1993). Elastic solids with many cracks and related problems. *Advances in Applied Mechanics*, 30, 259–445. https://doi.org/10.1016/S0065-2156(08)70176-5

Karabulut, H., & Bouchon, M. (2007). Spatial variability and non-linearity of strong ground motion near a fault. *Geophysical Journal International*, 170(1), 262–274. https://doi.org/10.1111/j.1365-246X.2007.03406.x

Karrech, A., Regenauer-Lieb, K., & Poulet, T. (2011). A damaged visco-plasticity model for pressure and temperature sensitive geomaterials. *International Journal of Engineering Science*, 49(10), 1141–1150. https://doi.org/10.1016/j.ijengsci.2011.05.005

Kaus, B. J. (2010). Factors that control the angle of shear bands in geodynamic numerical models of brittle deformation. *Tectonophysics*, 484(1–4), 36–47. https://doi.org/10.1016/j.tecto.2009.08.042

Kranz, R. L. (1979). Crack growth and development during creep of Barre granite. *International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts*, 16(1), 23–35. https://doi.org/10.1016/0148-9062(79)90772-1

Lavier, L. L., & Buck, W. R. (2002). Half graben versus large-offset low-angle normal fault: Importance of keeping cool during normal faulting. *Journal of Geophysical Research*, 107(B6), ETG–8. https://doi.org/10.1029/2001jb000513

Lavier, L. L., Buck, W. R., & Poliakov, A. N. B. (2000). Factors controlling normal fault offset in an ideal brittle layer. *Journal of Geophysical Research*, 105(B10), 23431–23442. https://doi.org/10.1029/2000JB00108

Lemiale, V., Mühlhaus, H.-B., Moresi, L., & Stafford, J. (2008). Shear banding analysis of plastic models formulated for incompressible viscous flows. *Physics of the Earth and Planetary Interiors*, 171(1–4), 177–186. https://doi.org/10.1016/j.pepi.2008.07.038

Le Pourhiet, L. (2013). Strain localization due to structural softening during pressure sensitive rate independent yielding. *Bulletin de la Societe Geologique de France*, 184(4–5), 357–371. https://doi.org/10.2113/gsgfbull.184.4-5.357

Lockner, D. A. (1998). A generalized law for brittle deformation of Westerly granite. *Journal of Geophysical Research*, 103(B3), 5107–5123. https://doi.org/10.1029/97JB03211

Lockner, D. A., Byerlee, J. D., Kuksenko, V., Ponomarev, A., & Sidorin, A. (1991). Quasi-static fault growth and shear fracture energy in granite. *Nature*, 350(6313), 39–42. https://doi.org/10.1038/350039a0

Lyakhovsky, V., Ben-Zion, Y., & Agnon, A. (1997). Distributed damage, faulting, and friction. *Journal of Geophysical Research*, 102(B12), 27635–27649. https://doi.org/10.1029/97JB01896

Manaker, D. M., Turcotte, D. L., & Kellogg, L. H. (2006). Flexure with damage. *Geophysical Journal International*, 166(3), 1368–1383. https://doi.org/10.1111/j.1365-246X.2006.03067.x

Marrett, R., Ortega, O. J., & Kelsey, C. M. (1999). Extent of power-law scaling for natural fractures in rock. *Geology*, 27(9), 799. https://doi.org/10.1130/0091-7613(1999)027(0799:EOPLSF)2.3.CO;2

McBeck, J., Kandula, N., Aiken, J. M., Cordonnier, B., & Renard, F. (2019). Isolating the factors that govern fracture development in rocks throughout dynamic in situ X-ray tomography experiments. *Geophysical Research Letters*, 46(20), 11127–11135. https://doi.org/10.1029/2019GL084613

Menéndez, B., Zhu, W., & Wong, T.-F. (1996). Micromechanics of brittle faulting and cataclastic flow in Berea sandstone. *Journal of Structural Geology*, 18(1), 1–16. https://doi.org/10.1016/0191-8141(95)00076-P

Meyer, S. E., Kaus, B. J. P., & Passchier, C. (2017). Development of branching brittle and ductile shear zones: A numerical study. *Geochemistry, Geophysics, Geosystems*, 18(6), 2054–2075. https://doi.org/10.1002/2016GC006793

Molnar, P., Anderson, R. S., & Anderson, S. P. (2007). Tectonics, fracturing of rock, and erosion. *Journal of Geophysical Research*, 112(F3). https://doi.org/10.1029/2005JF000433

Moresi, L., Dufour, F., & Mühlhaus, H.-B. (2003). A Lagrangian integration point finite element method for large deformation modeling of viscoelastic geomaterials. *Journal of Computational Physics*, 184(2), 476–497. https://doi.org/10.1016/S0021-9991(02)00031-1

Morley, C. K. (1996). Discussion of potential errors in fault heave methods for extension estimates in rifts, with particular reference to fractal fault populations and inherited fabrics. *Geological Society, London, Special Publications*, 99(1), 117–134. https://doi.org/10.1144/GSL.SP.1996.099.01.10

Naliboff, J. B., Glerum, A., Brune, S., Péron-Pinvidic, G., & Wrona, T. (2020). Development of 3-d rift heterogeneity through fault network evolution. *Geophysical Research Letters*, 47(13), e2019GL086611. https://doi.org/10.1029/2019GL086611

Nemat-Nasser, S., & Horii, H. (1982). Compression-induced nonplanar crack extension with application to splitting, exfoliation, and rockburst. *Journal of Geophysical Research*, 87(B8), 6805–6821. https://doi.org/10.1029/JB087iB08p06805

Olive, J.-A., Behn, M. D., Mittelstaedt, E., Ito, G., & Klein, B. Z. (2016). The role of elasticity in simulating long-term tectonic extension. *Geophysical Journal International*, 205(2), 728–743. https://doi.org/10.1093/gji/ggw044
Pan, S., Naliboff, J., Bell, R., & Jackson, C. (2022). Bridging spatiotemporal scales of normal fault growth during continental extension using high-resolution 3d numerical models. *Geochemistry, Geophysics, Geosystems*, 23(7), e2021GC010316. https://doi.org/10.1029/2021GC010316

Pan, S., Naliboff, J., Bell, R., & Jackson, C. (2023). How do rift-related fault network distributions evolve? Quantitative comparisons between natural fault observations and 3d numerical models of continental extension. *Tectonics*, 42(10), e2022TC007659. https://doi.org/10.1029/2022TC007659

Paterson, M. S., & Wong, T.-F. (2005). In *Experimental rock deformation: The brittle field* (Vol. 348). Springer.

Perol, T., & Bhat, H. S. (2016). Micromechanics-based permeability evolution in brittle materials at high strain rates. *Pure and Applied Geophysics*, 173(8), 2857–2868. https://doi.org/10.1007/s00024-016-1354-4

Petit, L., & Olive, J.-A. (2024). A sub-critically altered Maxwell (SCAM) constitutive law for brittle rocks: January 6, 2024 release (version 2.0) [Software]. Zenodo. https://doi.org/10.5281/zenodo.10464422

Poliakov, A. N. B., & Buck, W. R. (1998). Mechanics of stretching elastic-plastic-viscous layers: Applications to slow-spreading mid-ocean ridges. In W. Roger Buck, P. T. Delaney, J. A. Karson, & Y. Lagabrielle (Eds.), *Geophysical Monograph Series* (pp. 305–324). American Geophysical Union.

Rackauckas, C., & Nie, Q. (2017). DifferentialEquations.jl—A performant and feature-rich ecosystem for solving differential equations in Julia. *Journal of Open Research Software*, 5(1), 15. https://doi.org/10.5334/jors.151

Rudnicki, J., & Rice, J. (1975). Conditions for the localization of deformation in pressure-sensitive dilatant materials. *Journal of the Mechanics and Physics of Solids*, 23(6), 371–394. https://doi.org/10.1016/0022-5096(75)90001-0

Sammis, C. G., & Ashby, M. F. (1986). The failure of brittle porous solids under compressive stress states. *Acta Metallurgica*, 34(3), 511–526. https://doi.org/10.1016/0001-6160(86)90087-8

Scholz, C., Dawers, N., Yu, J.-Z., Anders, M., & Cowie, P. (1993). Fault growth and fault scaling laws: Preliminary results. *Journal of Geophysical Research*, 98(B12), 21951–21961. https://doi.org/10.1029/93jb01008

Spiegelman, M., May, D. A., & Wilson, C. R. (2016). On the solvability of incompressible Stokes with viscoplastic rheologies in geodynamics: Issues in viscoplasticity. *Geochemistry, Geophysics, Geosystems*, 17(6), 2213–2238. https://doi.org/10.1002/2015GC006228

Tada, H., Paris, P. C., & Irwin, G. R. (1973). The stress analysis of cracks. *Handbook, Del Research Corporation*, 34, 635.

Tapponnier, P., & Brace, W. F. (1976). Development of stress-induced microcracks in Westerly granite. *International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts*, 13(4), 103–112.

Tarantola, A. (2005). Inverse problem theory and methods for model parameter estimation. *Society for Industrial and Applied Mathematics*. https://doi.org/10.1137/1.9780898717921

Thomas, M. Y., Bhat, H. S., & Klinger, Y. (2017). Effect of brittle off-fault damage on earthquake rupture dynamics. In M. Y. Thomas, T. M. Mitchell, & H. S. Bhat (Eds.), *Geophysical Monograph Series* (pp. 255–280). John Wiley & Sons, Inc. https://doi.org/10.1002/9781119156895.ch14

Tsitouras, C., Simos, T. E., Psihoyios, G., & Tsitouras, C. (2009). Runge-Kutta pairs of orders 5(4) using the minimal set of simplifying assumptions. In *AIP Conference Proceedings* (pp. 69–72). Rethymno, Crete (Greece): AIP. https://doi.org/10.1063/1.3241561

Vermeer, P. A., & De Borst, R. (1984). Non-associated plasticity for soils, concrete and rock. *Heron*, 29(3).

Walsh, J. B. (1965a). The effect of cracks on the compressibility of rock. *Journal of Geophysical Research*, 70(2), 381–389. https://doi.org/10.1029/JZ070i002p00381

Walsh, J. B. (1965b). The effect of cracks on the uniaxial elastic compression of rocks. *Journal of Geophysical Research*, 70(2), 399–411. https://doi.org/10.1029/JZ070i002p00399

Wang, X.-Q., Schubnel, A., Fortin, J., Guéguen, Y., & Ge, H.-K. (2013). Physical properties and brittle strength of thermally cracked granite under confinement: Strength of thermally cracked granite. *Journal of Geophysical Research: Solid Earth*, 118(12), 6099–6112. https://doi.org/10.1002/2013JB010340

Wawersik, W. R., & Brace, W. F. (1971). Post-failure behavior of a granite and diabase. *Rock Mechanics Felsmechanik Mcanique des Roches*, 3(2), 61–85. https://doi.org/10.1007/BF01239627

Wu, C., Peng, Z., & Ben-Zion, Y. (2009). Non-linearity and temporal changes of fault zone site response associated with strong ground motion. *Geoophysical Journal International*, 176(1), 265–278. https://doi.org/10.1111/j.1365-246X.2008.04005.x

Wu, X., Baud, P., & Wong, T. (2000). Micromechanics of compressive failure and spatial evolution of anisotropic damage in Darley Dale sandstone. *International Journal of Rock Mechanics and Mining Sciences*, 37(1–2), 143–160. https://doi.org/10.1016/S1365-1609(99)00093-3

Zhang, J., Wong, T.-F., & Davis, D. M. (1990). Micromechanics of pressure-induced grain crushing in porous rocks. *Journal of Geophysical Research*, 95(B1), 341–352. https://doi.org/10.1029/JB095iB01p00341

Zoback, M. D., & Byerlee, J. D. (1975). The effect of microcrack dilatancy on the permeability of westerly granite. *Journal of Geophysical Research*, 80(5), 752–755. https://doi.org/10.1029/JB080i005p00752
Pan, S., Naliboff, J., Bell, R., & Jackson, C. (2022). Bridging spatiotemporal scales of normal fault growth during continental extension using high-resolution 3d numerical models. *Geochemistry, Geophysics, Geosystems*, 23(7), e2021GC010316. https://doi.org/10.1029/2021GC010316

Pan, S., Naliboff, J., Bell, R., & Jackson, C. (2023). How do rift-related fault network distributions evolve? Quantitative comparisons between natural fault observations and 3d numerical models of continental extension. *Tectonics*, 42(10), e2022TC007659. https://doi.org/10.1029/2022TC007659

Paterson, M. S., & Wong, T.-F. (2005). In *Experimental rock deformation: The brittle field* (Vol. 348). Springer.

Perol, T., & Bhat, H. S. (2016). Micromechanics-based permeability evolution in brittle materials at high strain rates. *Pure and Applied Geophysics*, 173(8), 2857–2868. https://doi.org/10.1007/s00024-016-1354-4

Petit, L., & Olive, J.-A. (2024). A sub-critically altered Maxwell (SCAM) constitutive law for brittle rocks: January 6, 2024 release (version 2.0) [Software]. Zenodo. https://doi.org/10.5281/zenodo.10464422

Poliakov, A. N. B., & Buck, W. R. (1998). Mechanics of stretching elastic-plastic-viscous layers: Applications to slow-spreading mid-ocean ridges. In W. Roger Buck, P. T. Delaney, J. A. Karson, & Y. Lagabrielle (Eds.), *Geophysical Monograph Series* (pp. 305–324). American Geophysical Union.

Rackauckas, C., & Nie, Q. (2017). DifferentialEquations.jl—A performant and feature-rich ecosystem for solving differential equations in Julia. *Journal of Open Research Software*, 5(1), 15. https://doi.org/10.5334/jors.151

Rudnicki, J., & Rice, J. (1975). Conditions for the localization of deformation in pressure-sensitive dilatant materials. *Journal of the Mechanics and Physics of Solids*, 23(6), 371–394. https://doi.org/10.1016/0022-5096(75)90001-0

Sammis, C. G., & Ashby, M. F. (1986). The failure of brittle porous solids under compressive stress states. *Acta Metallurgica*, 34(3), 511–526. https://doi.org/10.1016/0001-6160(86)90087-8

Scholz, C., Dawers, N., Yu, J.-Z., Anders, M., & Cowie, P. (1993). Fault growth and fault scaling laws: Preliminary results. *Journal of Geophysical Research*, 98(B12), 21951–21961. https://doi.org/10.1029/93jb01008

Spiegelman, M., May, D. A., & Wilson, C. R. (2016). On the solvability of incompressible Stokes with viscoplastic rheologies in geodynamics: Issues in viscoplasticity. *Geochemistry, Geophysics, Geosystems*, 17(6), 2213–2238. https://doi.org/10.1002/2015GC006228

Tada, H., Paris, P. C., & Irwin, G. R. (1973). The stress analysis of cracks. *Handbook, Del Research Corporation*, 34, 635.

Tapponnier, P., & Brace, W. F. (1976). Development of stress-induced microcracks in Westerly granite. *International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts*, 13(4), 103–112.

Tarantola, A. (2005). Inverse problem theory and methods for model parameter estimation. *Society for Industrial and Applied Mathematics*. https://doi.org/10.1137/1.9780898717921

Thomas, M. Y., Bhat, H. S., & Klinger, Y. (2017). Effect of brittle off-fault damage on earthquake rupture dynamics. In M. Y. Thomas, T. M. Mitchell, & H. S. Bhat (Eds.), *Geophysical Monograph Series* (pp. 255–280). John Wiley & Sons, Inc. https://doi.org/10.1002/9781119156895.ch14

Tsitouras, C., Simos, T. E., Psihoyios, G., & Tsitouras, C. (2009). Runge-Kutta pairs of orders 5(4) using the minimal set of simplifying assumptions. In *AIP Conference Proceedings* (pp. 69–72). Rethymno, Crete (Greece): AIP. https://doi.org/10.1063/1.3241561

Vermeer, P. A., & De Borst, R. (1984). Non-associated plasticity for soils, concrete and rock. *Heron*, 29(3).

Walsh, J. B. (1965a). The effect of cracks on the compressibility of rock. *Journal of Geophysical Research*, 70(2), 381–389. https://doi.org/10.1029/JZ070i002p00381

Walsh, J. B. (1965b). The effect of cracks on the uniaxial elastic compression of rocks. *Journal of Geophysical Research*, 70(2), 399–411. https://doi.org/10.1029/JZ070i002p00399

Wang, X.-Q., Schubnel, A., Fortin, J., Guéguen, Y., & Ge, H.-K. (2013). Physical properties and brittle strength of thermally cracked granite under confinement: Strength of thermally cracked granite. *Journal of Geophysical Research: Solid Earth*, 118(12), 6099–6112. https://doi.org/10.1002/2013JB010340

Wawersik, W. R., & Brace, W. F. (1971). Post-failure behavior of a granite and diabase. *Rock Mechanics Felsmechanik Mcanique des Roches*, 3(2), 61–85. https://doi.org/10.1007/BF01239627

Wu, C., Peng, Z., & Ben-Zion, Y. (2009). Non-linearity and temporal changes of fault zone site response associated with strong ground motion. *Geoophysical Journal International*, 176(1), 265–278. https://doi.org/10.1111/j.1365-246X.2008.04005.x

Wu, X., Baud, P., & Wong, T. (2000). Micromechanics of compressive failure and spatial evolution of anisotropic damage in Darley Dale sandstone. *International Journal of Rock Mechanics and Mining Sciences*, 37(1–2), 143–160. https://doi.org/10.1016/S1365-1609(99)00093-3

Zhang, J., Wong, T.-F., & Davis, D. M. (1990). Micromechanics of pressure-induced grain crushing in porous rocks. *Journal of Geophysical Research*, 95(B1), 341–352. https://doi.org/10.1029/JB095iB01p00341

Zoback, M. D., & Byerlee, J. D. (1975). The effect of microcrack dilatancy on the permeability of westerly granite. *Journal of Geophysical Research*, 80(5), 752–755. https://doi.org/10.1029/JB080i005p00752