Grain-size effects during semi-brittle flow of calcite rocks

We study the role of grain size in the rheological behaviour of calcite aggregates in the semi-brittle regime. We conduct triaxial deformation tests on three rocks, Solnhofen limestone, Carrara marble and Wombeyan marble, with average grain sizes of 5-10 $\mu$m, 200 $\mu$m and 2 mm, respectively, at pressures in the range 200-800 MPa and temperatures in the range 20-400 $^\circ$C. At all conditions, both strength and hardening rate increase with decreasing grain size. Flow stress scales with the inverse of grain size to a power between 1/3 and 2/3. Hardening rate decreases linearly with the logarithm of grain size. In-situ ultrasonic monitoring reveals that P-wave speed tends to decrease with increasing strain, and that this decrease is more marked at room temperature than at 200 and 400 $^\circ$C. The decrease in wave speed is consistent with microcracking, which is more prevalent at low temperature and low pressure. Microstructural observations reveal high twin densities in all deformed samples. Twin density increases with stress, consistent with previous datasets. Spatial distributions of intragranular misorientation indicate that twins are sometimes obstacles to dislocation motion, but this effect is not ubiquitous. Computed slip-transfer statistics indicate that that twins are typically weaker barriers to dislocation glide than grain boundaries, so that their effect on dislocation accumulation and hardening rates is likely smaller than the effect of grain size. Indeed, our data reveal that grain size exerts a first-order control on flow stress and hardening in calcite, whereas twinning may only have a secondary impact on these behaviours.


Introduction
In the shallow lithosphere, rocks deform by localised brittle failure and strength is described by a friction law (Scholz, 2002;Townend & Zoback, 2000). In the lower lithosphere, high temperatures and pressures promote the onset of crystal-plastic deformation mechanisms. Here, strength is described by flow laws sensitive to temperature (T ) and strain rate (ε) (Goetze & Brace, 1972;Evans & Kohlstedt, 1995). At intermediate conditions, deformation is ductile, i.e., remains macroscopically distributed (following the terminology of Rutter (1986)), and is often termed "semi-brittle", which is characterised by coupled cracking and crystal plasticity. Semi-brittle deformation is likely to support the highest stresses in the lithosphere (Goetze & Evans, 1979;Brace & Kohlstedt, 1980), impacting geodynamic processes (Burov, 2011) and the maximum depth of earthquake nucleation (Scholz, 1998). Despite the significance of semi-brittle flow, there is a paucity of simple models to describe the rheological behaviour of rock in this regime, and the relative importance of key processes, such as friction, tensile cracking, dislocation motion, twinning, or grain-boundary sliding is not well constrained. Thus, to advance our understanding of semi-brittle flow, we must first improve our understanding of interactions among microscale deformation processes.
Calcite rock is an important constituent of tectonic terrains, and undergoes a transition to semi-brittle flow at modest pressure and low temperature (e.g., Heard, 1960;Fredrich et al., 1989). This combination of characteristics has led to numerous laboratory investigations into the rheological behaviour of calcite rock across the brittle-ductile transition (see Rybacki et al. (2021) and references therein). The main deformation behaviours of calcite rock can be illustrated using the case of Carrara marble, a relatively isotropic, pure calcite aggregate with little initial crack porosity and equant grain shapes with sizes in the range 60-200 µm. At room temperature, Carrara marble is dilatant and brittle at low pressure (P < 30 MPa), and its strength is controlled by friction. At intermediate pressures (30 < P < 300 MPa), calcite rocks deform in a semi-brittle manner, and strength becomes decreasingly pressure sensitive with increasing pressure (Fredrich -3-manuscript submitted to JGR: Solid Earth et al., 1989). When pressure is high enough (P > 300 MPa), strength becomes independent of confining pressure and deformation is nondilatant (Edmond & Paterson, 1972). Furthermore, increases in temperature reduce the pressure required to promote semi-brittle flow (Rybacki et al., 2021). Edmond and Paterson (1972) and Fischer and Paterson (1989) report in-situ changes of volumetric strain during deformation of calcite rocks, finding that dilatancy is reduced by increasing pressure. Furthermore, sample volume remains constant when strength is independent of pressure (Edmond & Paterson, 1972). More recently, Schubnel et al. (2005) and Rybacki et al. (2021) measured in-situ P-wave velocity (V p ) and demonstrated that wave-speed decreases during deformation are suppressed by pressure increases. In-situ measurements therefore support microstructural observations in implying that increasing pressure suppresses cracking and frictional sliding whilst promoting crystal-plastic processes.
Microstructural observations from calcite rocks deformed in the semi-brittle regime reveal that distributed cracking, twinning, and dislocation motion act together to accommodate strain (Fredrich et al., 1989). Post-mortem crack density decreases with increasing confining pressure (Fredrich et al., 1989), approaching zero when strength is independent of confining pressure. Additionally, deformation twins are present in calcite rocks deformed at conditions below 800 • C (Rybacki et al., 2021), and twin spacing depends on stress (Rowe & Rutter, 1990). Detailed strain measurements at the grain scale reveal that twinning occurs readily in well oriented grains (those with high Schmid factor for twinning), and that it is likely associated with a local backstress that causes hardening of twinned grains (Spiers, 1979). In addition, microscale strain mapping also indicates the existence of shear localised along grain boundaries at temperatures < 800 • C, potentially identifying grain-boundary sliding as a possible deformation mechanism (Quintanilla-Terminel et al., 2017).
Another common rheological behaviour of semi-brittle flow in calcite rocks is strain hardening, which can persist to high temperatures (< 800 • C, see Rybacki et al., 2021).
At low temperature, strain hardening can arise from microscopic frictional slip across distributed defects, such as grain boundaries (e.g., David et al., 2020). In dislocationmediated deformation regimes, strain hardening occurs due to an increase in dislocation density caused by inefficient recovery mechanisms (e.g., lack of dislocation climb) (Mecking & Kocks, 1981) and is sensitive to microstructure. In particular, grain boundaries can act as barriers to dislocation motion, resulting in increases in strength and sometimes -4-manuscript submitted to JGR: Solid Earth increases in hardening rates with decreasing grain size for most metals (grain-size strengthening, see Cordero et al., 2016), and some geological materials (e.g., olivine, . Alternatively, strain hardening in calcite rocks has been proposed to be controlled by twinning (Rybacki et al., 2021). Twins may not contribute significantly to the total strain, but they could indirectly control strength and hardening by providing additional barriers to dislocations. Indeed, TEM observations show that dislocation densities are elevated adjacent to twin boundaries (Barber & Wenk, 1979;Fredrich et al., 1989;Rybacki et al., 2021). In the metallurgy literature (De Cooman et al., 2018), the additional hardening provided by twins is commonly known as twinning-induced plasticity, or TWIP. TWIP originates from observations of high-manganese steels, where high hardening rates and ductility occur as a result of the formation of deformation twins (De Cooman et al., 2018). The hardening rates in these metals are attributed to a "dynamic Hall-Petch effect" in which twinning refines the intracrystalline microstructure and thereby reduces the dislocation mean free path. In this case, twin spacing largely controls strength and hardening and thereby limits sensitivity to grain size.
The semi-brittle regime in calcite rock is thus characterised by many interacting deformation mechanisms. One way to quantify the relative contribution of each mechanism to the overall rheological behaviour is to test the impact of independent variables beyond the usual pressure and temperature conditions and imposed strain rate. One such variable is the grain size. At room temperature, yield stress (σ y ) and the transition pressure to semi-brittle deformation increase with decreasing grain size (Olsson, 1974;Fredrich et al., 1990). Grain size also impacts the strength of marble tested at elevated temperature in the dislocation-creep regime, which could be consistent with a Hall-Petch effect (Renner et al., 2002). In the high-pressure, moderate-temperature range (P > 200 MPa, where twinning is ubiquitous, models of twinning-induced plasticity suggest that twin density rather than grain size is the main control on strength (Rybacki et al., 2021). To test this hypothesis, here we explore the role of grain size in semi-brittle rheological behaviour and the brittle-plastic transition in calcite rocks for comparison to the role of twin density.
We performed a series of experiments at a range of pressures (P ≤ 800 MPa) and temperatures (≤ 400 • C) using calcite rocks spanning two orders of magnitude in grain size: Solnhofen limestone (grain size less than 10 µm), Carrara marble (grain size on the order of 100 µm), and Wombeyan marble (grain size on the order of 1 mm). Experiments -5-manuscript submitted to JGR: Solid Earth Initial V p (m s −1 ) 5600 5900 are supplemented by in-situ measurements of axis-parallel P-wave speed and post-mortem microstructural investigations to infer deformation mechanisms. We find that grain size has a first-order impact on both strength and hardening rate, and that the influence of twins is only indirect and likely smaller than anticipated from TWIP models.

Sample materials
Three calcite rocks were selected for experiments, with grain size, D, spanning three orders of magnitude. Solnhofen limestone is a lithographic limestone with a grain size of 5-10 µm (French et al., 2022), composed of > 99.9% calcite, and with an initial porosity (measured with helium pycnometry) of 4% (Baud et al., 2000). Carrara Marble is a medium-grained marble with equant grains of size 60-220 µm, comprised of > 99.9% calcite with a porosity of < 0.1%. In the starting material, most grains exhibit at least one twin set. Wombeyan marble is a coarse-grained marble with a grain size of 1-2 mm, comprised of 96 % calcite. Grains are equant and typically twinned in at least one plane, with an initial twin density of 60 mm −1 . This material was obtained thanks to Ian Jackson at the Australian National University, and is the same rock that was used by Paterson (1958). See Table 1 for further petrographic details.

Mechanical testing
A total of 42 experiments were conducted in the "Murrell" gas-medium triaxial ap-  mm in diameter and 22 mm in length, were dried in an oven at 70 • C for at least 24 hours before testing. Samples were inserted into annealed copper jackets of 0.1 mm wall thickness and were swaged onto the deformation rams using a slip ring (Figure 1a). The jacketed samples were inserted into the pressure vessel and pressurised to the target confining pressure using argon gas. Here the confining pressure is provides the intermediate and minimum principal stresses. Throughout the tests, the measured pressure remained within 5% of the target pressure. An internal furnace was used to heat samples to 200 • C and 400 • C. Axial load was applied vertically, generating the differential stress, σ. Axial stress was measured by an external load cell, and axial shortening was measured by a pair of linear variable-displacement transducers. Axial deformation was then applied by a piston that moved at a constant imposed shortening rate of 0.22 µm s −1 , equivalent to a strain rate of 10 −5 s −1 . All samples were deformed to a total strain of 7.5%.
A detailed list of test conditions is given in Table 3.

Data processing
Several corrections were made to the mechanical data. Seal friction and jacket strength were subtracted from load measurements (as reported in Harbord et al. (2022)) and displacement due to machine compliance was subtracted from shortening measurements.
The differential stress supported by the sample was computed by dividing the corrected force by the cross-sectional area of the sample, which was assumed to linearly increase with deformation consistent with a constant sample volume. Mechanical data were then further processed to obtain estimates of the tangent modulus (h = δσ/δε) by numerical differentiation of the stress data over a moving window of 1% strain. Uncertainties in seal friction during deformation result in error on the hardening modulus that is < 20 % of reported values. The yield stress (σ y ) is defined as the stress at which h falls to 90% of the sample-specific Young's modulus determined from elastic loading.

In-situ P-wave speed measurements
In-situ P-wave speed was measured during tests conducted on Solnhofen limestone and Carrara marble. Measurements were made parallel to the sample axis during deformation using the pulse-transmission method (Birch, 1961). Every 10 s during the experiments, a 200 V pulse of 0.4 µs duration and 2.5 MHz frequency was sent to a lead-titanatezirconate ceramic disk mounted centrally in the bottom ram, which was received by a transducer mounted centrally at the top of the sample assembly and recorded digitally at 100 MHz (see Harbord et al., 2022, for further details). The signal-to-noise ratio was improved by stacking 256 raw traces at each time interval. Changes in axial P-wave speed were computed relative to a reference waveform using cross correlation, and corrected for interfacial delays following the methods outlined in Harbord et al. (2022). Measurements of P-wave speed are reported as the change of wave speed (∆V ) and are normalised relative to the wave speed measured at the start of loading (V 0 ). Selected deformation tests were subsequently repeated, and wave-speed changes were found to be reproducible (see Supplementary Material Figure S1).
To gain insight into the microstructural state of samples after deformation during the return to ambient pressure, we continued to perform wave-speed surveys after removal of differential stress and throughout the staged decompression. The change in wave speed during decompression is quantified as the relative change in wave speed (∆V ) normalised by the wave speed measured at the end of deformation (V final ), and when at temperature, after cooling of the sample. These measurements were also complemented by postmortem measurements of wave speed at room pressure (0.1 MPa).

Optical microscopy
To investigate deformation microstructures, selected samples were mounted in epoxy, sectioned parallel to the deformation axis, and polished. A set of thin sections was also made for visible-light microscope observations. Imaging using transmitted and reflected light was performed using a Leica DM750P microscope furnished with an ICC50 camera. To quantify the prevalence of intragranular cracks in each sample, we followed the method outlined by Fredrich et al. (1989). Samples were imaged using reflected light at ×4 magnification, and we counted intersections between cracks (excluding grain boundaries) and a square grid of 2 mm by 2 mm with a spacing of 0.2 mm. These counts were used to determine the resulting crack surface area per volume, S v (mm 2 /mm 3 ).

Electron microscopy
Scanning electron microscopy was performed using a Jeol JSM-6480LV scanning electron microscope (SEM) hosted at University College London and a Zeiss Gemini 300 -9-manuscript submitted to JGR: Solid Earth field-emission gun SEM at the University of Cambridge. Electron backscatter diffraction (EBSD) patterns and forescattered electron images were collected using an Oxford Instruments Symmetry detector and AZtec 4.0 acquisition software. Table 2 lists the acquisition parameters for each EBSD dataset. All diffraction patterns with acquired with low detector gain. Datasets collected for conventional EBSD were acquired with a reduced number of pixels in the diffraction patterns to increase the mapping speed. Datasets collected for high-angular resolution electron backscatter diffraction (HR-EBSD) were collected with the maximum number of pixels permitted by the detector.
HR-EBSD is a postprocessing technique that analyses distortion of the diffraction patterns to measure lattice rotations and intragranular heterogeneity in elastic strain (Wilkinson et al., 2006;Britton & Wilkinson, 2011, 2012. Full details of the technique are given by Wallis et al. (2019) and here we provide a summary of the key points. One diffraction pattern from the host grain within each mapped area was manually selected to be a reference pattern based on the quality of the diffraction pattern and its position within the map. 100 regions of interest, 256×256 pixels in size, were extracted from all diffraction patterns within the host grain. Each region of interest from each diffraction pattern was cross-correlated with the corresponding region of interest from the reference pattern to determine shifts in their positions. Shifts in the diffraction pattern due to beam scanning were corrected using a calibration determined on an undeformed Si single crystal following Wilkinson et al. (2006) and the position of the pattern centre was calibrated using diffraction patterns collected over a range of detector insertion distances (Maurice et al., 2011). A deformation-gradient tensor was fit to the field of shifts in each diffraction pattern. The deformation-gradient tensor was decomposed into its symmetric and antisymmetric parts, which respectively give the elastic-strain and rotation tensors (Wilkinson et al., 2006). We used the pattern remapping approach of Britton and Wilkinson (2012), in which a first pass of cross-correlation measures lattice rotations that are used to rotate each pattern back into the orientation of the reference pattern before a second pass of cross-correlation measures the elastic strain and a small correction to the rotations.
The elastic strains were converted to stresses using Hooke's law. The measured stresses are relative to the unknown stress state at the reference point, giving maps of intragranular stress heterogeneity, rather than absolute values. We subtracted the arithmetic mean value of each component of the stress tensor within the map area from each measured value so that the final maps provide stress heterogeneity relative to the unknown mean -10-manuscript submitted to JGR: Solid Earth stress state within each grain (Mikami et al., 2015). Alongside, spatial gradients in the lattice rotations were used to estimate densities of geometrically necessary dislocations (GNDs). Densities of each dislocation type on the slip systems summarised by J. H. De Bresser and Spiers (1997) were fit to the measurable components of the lattice curvature following the approach applied to quartz by Wallis et al. (2019). We emphasise that the stress heterogeneity and GND densities are determined independently, being respectively derived from the distinct symmetric and antisymmetric parts of the deformationgradient tensor. Data points were filtered out if they had a mean normalised peak height in the cross-correlation function of <0.3 or a mean angular error in the fitted deformation gradient tensor of >0.004 radians (Britton & Wilkinson, 2011).
We analysed the probability distributions of the stress heterogeneity to assess whether the stresses are imparted by dislocations. The probability distribution of the stress field of population of dislocations has a characteristic form with tails that depart from a normal distribution towards higher stresses (Jiang et al., 2013;Wilkinson et al., 2014). We assess for the presence of these tails using a normal-probability plot, in which the cumulativeprobability axis is scaled such that a normal distribution falls on a straight line. Importantly, the tails of the probability distribution, P (σ), have a specific form if the stresses, σ, are imparted by dislocations, whereby P (σ) → Cρ|σ| −3 , where C is a constant that depends on the material, type(s) of dislocation, and considered stress component, and ρ is the total dislocation density (Groma & Bakó, 1998;Wilkinson et al., 2014). To test whether the measured stress fields exhibit this form, we compute the restricted second moment, ν 2 , which is a metric that characterises the shape of a probability distribution based on the integral over restricted ranges in stress, calculated as ν 2 (σ) = +σ −σ P (σ)σ 2 dσ (Wilkinson et al., 2014;Kalácska et al., 2017). A plot of ν 2 versus ln(σ) becomes a straight line at high stresses if the tails of the probability distribution of the stresses exhibit the form P (σ) ∝ |σ| −3 expected of a population of dislocations (Wilkinson et al., 2014;Kalácska et al., 2017). We apply this analysis to the σ 12 component of the stress tensor as this component is the least modified by sectioning the sample and is a shear stress capable of exerting glide forces on dislocations . This approach has recently been applied to olivine by Wallis et al. (2021Wallis et al. ( , 2022. We include in these plots data from an undeformed Si wafer measured by Wallis et al. (2022) to provide an indication of the noise level of the stress measurements.

Twin density measurements
Twin density was measured using a combination of forescattered electron images and EBSD maps. For each sample for which twin density was reported, we chose between 20 and 60 grains in a representative forescattered electron image, and measured twin spacing and twin width perpendicular to the selected twin set. Using an EBSD map of the same area, we determined the orientation of each grain and active twin set, and used the angle between the normal to the twin plane and the normal to the section to correct the measured twin width and spacing (the same procedure as used by Rutter et al., 2022).

Effect of pressure and temperature
The absolute strength, i.e. the differential stress at a given strain, of all rocks tested decreases with increasing temperature. The temperature sensitivity is greater at elevated pressure (400 and 600 MPa) than at low pressure, and is greater in Wombeyan marble than in Carrara marble and Solnhofen limestone ( Figure 3).
Both Carrara marble and Wombeyan marble have strengths and yield stresses that tend to increase with increasing pressure ( Figure 3). However, strength seems to become pressure-independent beyond a pressure threshold. In Carrara marble, there are no quantitative differences between the tests conducted at 600 MPa and at 800 MPa, and only small differences can be detected between 400 MPa and 600 MPa. In Wombeyan marble, strength does not change appreciably with pressure above 400 MPa. The pressure sensitivity of strength in Carrara marble and Wombeyan marble is reduced by increasing temperature, and strength becomes pressure independent at a lower pressure as temperature increases. Hardening rates in Carrara marble and Wombeyan marble also increase with increasing pressure, except for at a temperature of 400 • C, at which hardening rates are independent of pressure.
Solnhofen limestone exhibits some differences compared to the other lithologies in terms of the pressure dependence of strength. At all conditions, the yield stress of Solnhofen limestone decreases with increasing pressure, in contrast with the other lithologies. In  Figure 4. Evolution of yield stress (σy), differential stress at 2.5% strain (σ2.5) and 5% strain (σ5) as functions of grain size. Circles correspond to data obtained in this work, and other symbols are taken from the literature (see Table 4 for list of labels).  Table 4 for list of symbols).
contrast, at 5% strain, the strength of Solnhofen limestone increases slightly with pressure from 200 to 400 MPa, but decreases at 600 MPa ( Figure 3a).
Wave-speed changes also depend on pressure and temperature. Decreases in axial P-wave speed with strain are smaller at higher temperature in both Solnhofen lime-

Effect of grain size
Grain-size has a marked effect on the mechanical behaviour. Solnhofen limestone is always the strongest material and Wombeyan marble the weakest, independent of test conditions (Figures 2 and 4). For example, at a temperature of 20 • C and pressure of 600 MPa, the stress at 7.5% strain is 570 MPa for Solnhofen limestone, 400 MPa for Carrara marble, and 310 MPa for Wombeyan marble (Figure 4).
Grain-size strengthening is further promoted by increases in strain. For example, at a temperature of 20 • C the range of stresses across the grain sizes tested is σ y = 100-300 MPa at the yield stress, σ 2.5 = 250-500 MPa at 2.5% strain, and σ 5 = 300-600 MPa at 5% -17-manuscript submitted to JGR: Solid Earth Table 4. Summary of references for literature data presented in Figure 4 and 5.
Wave-speed evolution is also sensitive to grain size, as highlighted by comparing the behaviour of Solnhofen limestone and Carrara marble. The overall decrease in wave speed is nearly always greater for Solnhofen limestone at a given set of conditions than it is for Carrara marble. For example, at T = 20 • C, P = 600 MPa and ε = 7.5%, the wave-speed change is −5% in Carrara marble ( Figure 2b) but is −11% in Solnhofen limestone ( Figure 2a). The only exception to this grain-size dependence is at T = 400 • C and P = 600 MPa. At these conditions, the final relative wave-speed change at ε = 7.5% is +2% in Solnhofen limestone ( Figure 2g) and -2% in Carrara marble (Figure 2 h) dashed dark blue curve).

Wave-speed changes during decompression
Wave speed always decreases significantly during decompression ( Figure 6). In Solnhofen limestone, the P-wave speed remains approximately constant during initial decompression down to around 400 MPa. Further decompression leads to a decrease in wave speed, which drops substantially at pressures below 200 MPa. In samples deformed at room temperature, the wave speed of the recovered material is as low as 50% of the wave speed measured after deformation but prior to decompression. This drop becomes less marked in samples deformed at high temperature, with changes on the order of 25-35% at T = 200 • C, and 15-30% at T = 400 • C. The behaviour is generally similar in Carrara marble, with some notable quantitative differences. The wave speed starts decreasing immediately as pressure in decreased, even in tests conducted at P = 600 MPa ( Figure   6b,d). The characteristic pressure below which P-wave speed drops most markedly is on the order of 100 MPa. The effect of deformation temperature is similar to that observed in Solnhofen limestone, with elevated temperature during deformation promoting a more limited reduction in wave speed during decompression.

Brittle features
Systematic observations of cracks within deformed samples reveals three distinct deformation regimes according to deformation conditions (figure 7). The low-pressure, low-temperature regime and the high-pressure, low-temperature regime are characterised by cracks, whereas in the high-temperature regime cracks are absent.
In the low temperature regimes where cracks are observed, cracks interact with twins.
Dense arrays of cracks confined to individual twin lamellae are commonly observed (Fig--20-manuscript submitted to JGR: Solid Earth  -22-manuscript submitted to JGR: Solid Earth ure 8c). In addition, some intragranular cracks contain steps and change orientation across individual twin lamellae ('stair-step', Figure 8e). In other instances, microcracks are observed to nucleate from the tips of twins, and do not reach grain boundaries (Figure 8f).
Low-pressure, low-temperature regime The low-pressure, low-temperature regime is limited to Carrara marble experiments conducted at 200 MPa and 20 and 200 • C only. In this regime mode-I intragranular microcracks are a common feature, and result in the highest crack densities. These microcracks are typically confined to single grains and are of low aspect ratio (e.g., Figure 8b,d).
In places, the microcracks are concentrated in discrete regions as crush zones (Figure 8a).
In these regions, cracks often span a few grains and result in a locally elevated crack density.
Cracks are also observed to relate to the geometry of grains and grain boundaries.
Some cracks nucleate at geometric irregularities along grain boundaries. For example, steps in grain boundaries are often associated with short tensile cracks that propagate a short distance into the grain interior (Figure 8d). Smaller grains can also act as indenter grains, and cracks are nucleated to accommodate the indenter grain shape ( Figure   8b). Often, cracks relating to geometric incompatibilities can be wholly contained within a single grain (Figure 8b and d).
High-pressure, low-temperature regime The high-pressure, low-temperature regime spans samples deformed above P = Within the high-pressure, low-temperature regime, the crack density decreases with increasing pressure for Carrara marble and Wombeyan marble. For example in Carrara marble at 20 • C, S v = 7.31 mm 2 /mm 3 at P = 400 MPa reducing to S v = 2.28 mm 2 /mm 3 at P = 600 MPa (Table 5, Run0078 and Run0084). Furthermore, for a given set of con--23-manuscript submitted to JGR: Solid Earth ditions in the high-pressure, low-temperature regime, crack density is lower in the coarser grained Wombeyan marble.

Crystal-plastic microstructures
Crystal-plastic microstructural features are dominated by deformation twins that are present in samples deformed at all conditions ( Figure 9). Nearly all grains are twinned on at least one plane, grains that contain two twin sets are also common, and occasionally grains contain three twin sets. In particular, the density of twins, i.e., the number of lamellae per unit length, is high in samples deformed at low temperature (Figure 9 a-d) and decreases with increasing temperature (Figure 9g On a larger scale than twin lamellae, undulose extinction is widespread (Figure 9).
In Wombeyan marble, the intensity of the undulose extinction increases when approaching grain boundaries and grain-boundary irregularities (Figure 9c). In Wombeyan marble, the undulose extinction observed in some instances suggests corrugation of the crystal lattice (Figure 9e). In Carrara marble, grains also display undulose extinction, although it is not as common as in Wombeyan marble. Again, in Carrara marble undulose extinction is associated to geometric irregularities (Figure 9h).

SEM observations
Forescattered electron images were used to image twins at high spatial resolution.
These images reveal very thin deformation twins, especially in Solnhofen limestone, which in some places are on the order of 100 nm in thickness (Figure 10a,b). Grains of calcite in Solnhofen limestone also exhibit lower twin incidence than the coarser-grained samples of Wombeyan marble and Carrara marble (Figure 10c ( Figure 12c). In other grains, the GROD pattern is largely unaffected by twins and instead exhibits lattice distortion over larger length scales approaching the grain size (Figure 12d). In some cases, there is a mixed interaction, where the GROD weakly follows the twin orientation but is also affected by the grain geometry ( Figure 12b).
HR-EBSD maps ( Figure 13) reveal the distributions of GNDs and intragranular stress heterogeneity in the endmember grains in Figure 12. The grain in Figure 12b, which exhibited the least impact of the twins on the distribution of GROD, also exhibits negligible impact of the twins on GND density (Figure 13a). Within the grain interior, GND density is generally at or below the background noise level of approximately 10 13.5 m −2 arising from noise in the rotation measurements. Apparent GND densities above this noise level are highly localised in the vicinity of twin boundaries and potentially result from In addition to the EBSD mapping, based on the forescatter imaging, we were able to measure twin spacing and correct their thickness in the same manner as Rutter et al.
(2022, Figure 14). We report values of twin density, ρ t (mm −1 ), which is computed as the inverse of the measured average twin spacing (Figure 14b). Our data agree with previously reported measurements of stress versus twin density. Twin density is lowest in

Discussion
The mechanical data reveal that deformation of calcite rocks across the range of conditions that we tested is ductile. Almost all experiments demonstrate strain-hardening behaviour but the precise characteristics of the hardening vary with pressure and temperature. The deformation of calcite rocks in this study is accommodated by cracking, twinning and dislocation motion, depending on experimental conditions. The main question we seek to answer is what controls the rheological behaviour of calcite rocks in the semi-brittle regime? Determination of the rheological behaviour of calcite rocks requires knowledge of the approximate partitioning of strain between each active deformation mechanism. This can be answered by considering the microstructural, wave-speed and mechanical data.

Brittle deformation
The creasing hardening rates, smaller wave-speed decreases and a reduction in post-mortem crack density ( Figure 15 and Table 5). Similar observations of decreasing post-mortem crack density with increasing pressure were made by Fredrich et al. (1989) in Carrara marble at room temperature. Measurements of volumetric strain during the deformation of Carrara marble and Solnhofen limestone at room temperature made by Edmond and Paterson (1972) demonstrate that dilatancy is suppressed at high pressure. Both of these studies also report increased hardening rates with increasing pressure. Decreasing strain accommodation by cracking is therefore systematically associated with increasing hardening rates ( Figure 15).
With increasing temperature, pressure-independent strength is reached at a lower pressure. In conjunction with this, deformation at elevated pressure is systematically related to smaller drops in wave speed and increased hardening rates. Crack density is also lower for a given pressure at higher temperature. Measurements of reduced pore volume during deformation of Carrara marble at high temperatures by Fischer and Paterson (1989) support this observation. Taken together, these observations identify that cracking is suppressed by temperature increases.
Another interesting observation is the large wave-speed decrease during decompression. This feature suggests that significant sample damage is accrued during decompression. Samples with the largest wave-speed decrease are those deformed at high-pressure, low-temperature conditions, which correspond to a low intragranular fracture density (ta-  Edmond and Paterson (1972), who reported increases in volumetric strain during decompression in deformed Solnhofen limestone and Carrara marble, with the magnitude of volumetric strain change increasing with the final stress level. This phenomenon is likely to result from deformation caused by strain incompatibility between grains (Ashby, 1970). During deformation at high pressure, individual grains deform both plastically and elastically to accommodate the imposed strain in the sample, which results in heterogeneous internal stresses. Upon removal of the confining pressure, these internal stresses are no longer in equilibrium with the applied external forces and are likely to generate interface cracks between grains to accommodate the strain incompatibilities accrued during deformation.
Based on experiments at room temperature and moderate pressure, Fredrich et al.
(1990) discussed grain-size strengthening in the model framework of Horii and Nemat-Nasser (1986). Horii and Nemat-Nasser (1986) solved a wing-crack model coupled to a 'plastic zone', and their results predicted that coarser-grained materials should be more 'brittle', i.e., the ratio of differential stress to confining pressure at the brittle-ductile transition should decrease with increasing grain size. The results of Fredrich et al. (1990) showed that this ratio was independent of grain size, which they explained by considering that the plastic yield stress scaled as D −1/2 . This hypothesis is supported by our observations. More generally, yield stress is observed to scale inversely with grain size in metals (Cordero et al., 2016) and also in olivine .
Grain-size strengthening is a common phenomenon in metals and has received significant attention in the metallurgy literature (see Cordero et al. (2016) and Y. Li et al.  (2016) for reviews). Strengthening with decreasing grain size is usually termed the Hall-Petch effect (Hall, 1951;Petch, 1953), and is typically of the form where σ 0 is the intrinsic resistance of a lattice to dislocation motion, K is a materialdependent Hall-Petch coefficient and m is a dimensionless exponent. In the original observations and theoretical arguments of Hall (1951) and Petch (1953), the exponent m = 0.5. Values of m = 0.2-1 have been reported in metals (Dunstan & Bushby, 2014;Cordero et al., 2016;Y. Li et al., 2016).
We fitted equation 1 to our data using a least-squares regression and setting m = 0.5 (Table 4.2), similar to Renner et al. (2002).  (Hall, 1951;Petch, 1953); (2) the grain-boundary ledge model, in which grain-and subgrain-boundary irregularities emit forest dislocations that act as obstacles (see J. C. Li, 1963); (3) the plastic-strain model, in which the rate of increase in dislocation density with plastic strain is inversely proportional to grain size (Conrad et al., 1967); (4) the elastic-anisotropy model, in which interactions among elastically anisotropic grains require the introduction of GNDs, and smaller grains have relatively larger strain gradients and greater GND densities (Meyers & Ashworth, 1982). All of these models, except that of Meyers and Ashworth (1982), arrive at an expression similar to Equation 1. The coefficient K always includes the shear modulus, G, and the Burgers vector, b, as well as other geometric constants, and m = 0.5. However, fitting exercises have shown that for a wide range of metals, m = 0.5 and may be better fit by m = −1 or by the relationship ln d/d (Dunstan & Bushby, 2014;Y. Li et al., 2016).
A grain-size dependence of the yield stress suggests that the plastic-strain model 3 is not appropriate for calcite as a finite amount of macroscopic plastic strain is required to generate the Hall-Petch effect (Conrad et al., 1967;Ashby, 1970). As pointed out by Hansen et al. (2019), models 1, 2 and 4 also require finite plastic strain, although this may be sufficiently localised for the macroscopic behaviour still appear elastic prior to macroscopic yielding (Maaß & Derlet, 2018). The temperature dependence of the yield stress suggests that the models must include short-range dislocation interactions, as longrange interactions are elastic and therefore largely temperature insensitive .
In model 1, the Hall-Petch effect is attributed to pile-ups of dislocations at grain boundaries, and their role in promoting yield of neighbouring grains (Hall, 1951;Petch, 1953). In this case, larger grains can support longer pile-ups, which generate more in-  Li, 1963). As materials with finer grain sizes have greater grain-boundary areas, there are more potential sites for generation of dislocations, resulting in a higher dislocation density and therefore stress required for macroscopic plastic deformation. In our samples it is difficult to assess the grain-boundary ledge model. Some evidence may be that undulose extinction is often controlled by the geometry of grain boundaries (Figure 9), indicating that dislocation activity is influenced by grain boundaries, although it does not identify whether grain boundaries act as dislocations sources.
The elastic-anistropy model of Meyers and Ashworth (1982) relies on incompatibilities in elastic strain between neighbouring grains. When stress is applied, GNDs form to accommodate the elastic mismatch between grains. Our microstructural observations and wave-speed data may be compatible with this model. Open grains boundaries (Figure 8 c, e and g) and large wave-speed decreases during decompression ( Figure 6) suggest significant relaxation of internal stress during decompression and, in turn, this internal stress may result from elastic mismatch among grains and associated GND for-mation. EBSD observations also reveal significant lattice curvature ( Figure 12) and increases in GND density in the vicinty of grain boundaries ( Figure 13). Furthermore, microscale observations of strain in deformed calcite rocks by Spiers (1979) and Quintanilla-Terminel and Evans (2016) reveal local heterogeneity in finite strain near to grain boundaries. In particular, Quintanilla-Terminel and Evans (2016) identify strain heterogeneity on a scale similar to the grain size. Wallis et al. (2018) observed an increase in misorientation in the vicinity of grain boundaries in naturally deformed calcite rocks, consistent with the presence of plastic strain gradients imparted in response to strain incompatibility among neighbouring grains (Meyers & Ashworth, 1982).
In summary, our observations suggest that either model 1, 2 or 3 maybe applicable to calcite. More systematic microstructural observations are required to discriminate among these models. As a generality, all models predict that the strength is sensitive to the mean free path of dislocations, which is controlled by grain size. Therefore, Hall-Petch models are closely related to Kocks-Mecking-Estrin (KME) single state-variable models of strength, which incorporate the role of dislocation mean free path (Kocks, 1966(Kocks, , 1976Mecking & Kocks, 1981). To explore the KME model we must consider the role of twins, which is discussed in the next section.

Is TWIP compatible with the deformation of calcite rocks?
The activity of mechanical twinning in calcite is closely related to the magnitude of differential stress (Jamison & Spang, 1976;Spiers, 1979;Rowe & Rutter, 1990). In experiments, Rowe and Rutter (1990) demonstrated that the spacing of twins decreases with increasing stress, independent of temperature and grain size. Twins are observed in all our deformed samples, and twin density increases with differential stress consistently with the observations of Rowe and Rutter (1990); Rybacki et al. (2013);Rutter et al. (2022). Rybacki et al. (2013) also demonstrated that stress was proportional to the square root of twin density at low stresses (< 250 MPa). However, at high stress, the value of stress saturates with respect to twin density and the square root dependence breaks down (Rowe & Rutter, 1990, and Figure 14). Transmission Electron Microscope observations of calcite deformation twins indicate that twins often interact with dislocations (Barber & Wenk, 1979;Fredrich et al., 1989;Rybacki et al., 2013). Given these observations, Rybacki et al. (2021) argued that twin spacing linearly decreases as dislocation density increases, as flow stress is proportional to the square root of dislocation -37-manuscript submitted to JGR: Solid Earth density (Taylor, 1934). These observations suggest that twin spacing either is a consequence of, or controls, the strength of calcite rocks. Rybacki et al. (2021) argued that twin spacing directly controls the strain hardening rate, and by extension the strength, by drawing analogy to twinning induced plasticity steels (TWIP, De Cooman et al., 2018). Models of TWIP originate from high manganese steels, which exhibit high hardening rates (approx. 3% of the shear modulus) with respect to other steels (approx. 0.05 % of the shear modulus), as a consequence of mechanical twinning. The mechanism of TWIP originates from the abrupt changes in crystallographic orientation at twin boundaries, which can act as barriers to dislocation motion. Progressive twinning leads to dynamic refinement of the microstructure (De Cooman et al., 2018), in which finer twin spacing reduces the mean free path of dislocations (λ) and causes a dynamic Hall-Petch effect.
In phenomological models of TWIP in the metallurgical literature, strain hardening is attributed to a dynamic Hall-Petch effect resulting from progressive twinning (Bouaziz et al., 2008). The model of Bouaziz et al. (2008) is formulated by first considering the Taylor equation that relates shear flow stress (τ ) to the total dislocation density (ρ), where τ 0 is the initial strength of a polycrystal, τ b is the back stress, which may arise from long-range dislocation interactions and stress fields around twins, α a constant close to unity, b the Burgers vector and µ the shear modulus. In their formulation, the final term represents isotropic hardening due to short-range dislocation interactions. In calcite, the Taylor equation has been demonstrated to apply to calcite rocks deformed at temperatures of 550-700 • C (J. H. P. De Bresser, 1996), although the relative contributions of kinematic hardening due to long-range dislocation interactions that generate back stress and isotropic hardening due to short-range dislocation interactions have not been separated.
To obtain the evolution of stress with strain, the Taylor relation is combined with a modified Kocks-Mecking-Estrin equation (Kocks, 1966(Kocks, , 1976Mecking & Kocks, 1981) to describe the change of ρ with strain (Bouaziz et al., 2008): where f is a rate-and temperature-dependent dynamic-recovery coefficient, k is a constant that characterises dislocation storage due to dislocation interactions, and D t is the twin spacing. In this model, changes in the dislocation mean free path are sensitive to the total dislocation density, grain size and twin spacing. An additional term, given by the product of the recovery factor and the dislocation density, f ρ, is subtracted to account for dynamic recovery processes. Rybacki et al. (2021)  of the shear modulus) are not unique to calcite rocks as olivine, which does not twin, exhibits hardening rates up to 5-10 % of the shear modulus when deformed by low-temperature plasticity Druiventak et al., 2011).
The TWIP model also suggests that the hardening rate should be dominated by twin spacing as D t is always at least 1-2 orders of magnitude smaller than the grain size D (Fig. 14). Taking the case of Carrara marble at 600 MPa, 200 • C with D = 100 µm and D t = 5 µm, we might expect strain hardening of 450αµ (neglecting recovery and forest hardening). At the same conditions for Wombeyan marble, D = 1 mm and D t = 7.5 µm, so that strain hardening should be of 360αµ. The ratio of hardening rates would therefore be 1.25. The actual hardening rates observed in our experiments are in a ratio of 2 (considering H 5 = 1.6 and 0.8 GPa for Carrara marble and Wombeyan marble respectively), which suggests that twins may have a smaller impact on hardening than anticipated from Equation 3.

Twins as potential barriers to dislocations
Our main observation is that of grain size strengthening (Figure 4), and it is difficult to assess the exact role of twins in the hardening behaviour from our microscopic data alone. The TWIP model is founded on the notion that twins produce further hardening by either retarding or stopping dislocation motion. In order to test the potential validity of the TWIP model in calcite, in this Section we assess the respective efficiency of grain boundaries and twins at impeding dislocations.
Microstructural observations suggest that the interaction between dislocations and twins boundaries varies between grains. Some grains contain lattice curvature that is clearly affected by twin boundaries (e.g., Figure 12c) whereas in other grains the lattice curvature appears to be affected by twin boundaries (e.g., Figure 9c and Figure 12d). The trapping of dislocations therefore appears to depend on grain orientation and related twin orientation.
We can assess the effectiveness of twin boundaries and grain boundaries as barriers to dislocations by considering slip-transmission coefficients. The simplest form of this analysis is purely geometric and considers only the slip-plane orientation and direction of the Burgers vector of the incoming and outgoing slip systems (Luster & Morris, 1995),   (Table 4) demonstrate that m is always greater than 0.4, and in some cases twins are transparent with a value of 1 (e.g., r + 2 to r + 3 across an e 1 twin). Given also the large number of available slip systems, this analysis suggests that twin boundaries may impede dislocation motion less than randomly oriented grain boundaries. 0.596 To further compare the effect of twin and grain boundaries on dislocation motion, we can also consider the effects of twin-boundary orientation. This analysis can be per--41-manuscript submitted to JGR: Solid Earth formed with a geometric criterion that quantifies the degree of misalignment between the line intersections of incoming and outgoing slip planes with the twin plane, l A and l B , respectively, and between their slip directions, d A and d B , respectively. The following scalar quantity is maximal when the two slip systems on either side of the boundary are aligned and slip can be transmitted easily across the boundary (Shen et al., 1986;Bayerschen et al., 2016) The line intersections l can be obtained from l = (n × n Γ )/|n A × n Γ |, where n Γ is the twin boundary plane normal.
We used theM criterion to compare the efficiency of slip transmission across twin boundaries to slip transmission across grain boundaries. A random fabric was generated using MTEX, and we computedM between random grain pairs by assuming that the activated slip system (either r ± or f ± slip) was that with highest Schmid factor in each grain. For each grain, we also computedM across a twin boundary hosted in the initial grain, the activated twin system was assumed to be that with the highest Schmid factor in each grain. These calculations suggest an average value for slip transmission ofM = 0.27 across grain boundaries, which is considerably smaller than the average value ofM = 0.54 for slip transfer across twin boundaries.
To further test this result, we computed the expected m values for our EBSD data across the activated twin systems within each grain ( Figure 17). The results demonstrate that grains with low values of m ( Figure 17, grain c) exhibit dislocation substructures that are affected by twinning, that is segmentation of lattice distortion between twins ( Figure 12c). In contrast, grains with large values of m ( Figure 17, grain d) exhibit gradients in lattice distortion over length scales approaching the grain size, but which are not strongly affected at smaller scales by twin boundaries. HR-EBSD maps taken from the interior of these grains supports this conclusion, large residual stresses and GND densities are observed in the vicinity of twin boundaries in grains for which m is small (Figure 13b and c). Lower residual stresses and GND densities are present in the vicinity of twin boundaries in grains for which m is low (Figure 13a and Figure 17b). The efficacy of twin boundaries as barriers to dislocation motion is therefore dependent on the local orientation of individual grains, but is on average weaker than a grain boundary taken at random. The relative ease with which dislocations can transmit across twin boundaries suggests that Equation 3 requires refinement. We suggest that weights should be applied to the contributions of twin boundaries and grain boundaries, such that equation 3 becomes, in which k D and k t are weights to account for the relative efficacy of grain boundaries and twin boundaries. We expect that k t k D . Further microstructural measurements, such as the evolution of twin density with strain, are required to determine the value of these weighting factors.

Towards a model of semi-brittle flow in calcite rocks
Our observations combined with previous results suggest several key characteristics that should be captured by a model of semi-brittle flow in calcite-rich rocks at T < 400 • C: (1) non-linearly increasing strength and hardening with increasing pressure, (2) decreasing strength with increasing temperature, (3) increasing strength with decreasing grain size, (4) strength that is nearly insensitive to strain rate. Rybacki et al. (2021) reviewed proposed models of semi-brittle flow, identifying difficulties in combining brittle and plastic flow into a simple unified model.
Through the semi-brittle flow regime, the strain contribution of brittle and crystalplastic processes changes with depth. At low pressures and temperatures, at which brittle processes dominate, the macroscopic behaviour of rocks is described by frictional failure (Byerlee & Brace, 1966;Brace & Kohlstedt, 1980). Microphysical models of brittle deformation are typically based the wing-crack model (e.g., Nemat-Nasser & Horii, 1982;Ashby & Sammis, 1990), in which brittle damage is accounted for by the propagation of Mode I wing cracks. In this regime, strength is pressure sensitive, dependent on grain size and is to first order strain rate-insensitive. There is no strong temperature dependence, and plasticity is not considered.
There are only a small number of microphysical models in existence accounting for coupled brittle-plastic deformation. Horii and Nemat-Nasser (1986) modified and solved the problem of a wing crack coupled a 'plastic zone' by considering a dislocation pileup ahead of the shear crack. Their model is sensitive to pressure, grain size and also temperature as the plastic yield strength in the plastic zone can vary with temperature. How-ever, their model predicts that materials with coarser grain size are more brittle, contradicting observations from experimental data (Fredrich et al., 1990).
More recently, Nicolas et al. (2017) derived a model incorporating the propagation of wing cracks, plastic pore collapse and nucleation of new cracks due to dislocation pileups. The model reproduces several important characteristics of the deformation of porous limestones. As pointed out by Rybacki et al. (2021), this model does not consider dynamic recovery or twinning, which are important deformation processes in calcite-rich rocks. The large number of free parameters make this model challenging to test and may limit its general application.
An alternative approach to introduce feedbacks between cracking and plastic flow may be by use of a modified Kocks-Mecking equation. The anticorrelation between strain hardening and crack density indicated by our results (Figure 15) suggests that cracking acts to reduce stress and hardening. The role of cracks could be multiple. One possibility is that strength remains dictated by dislocation density, and correctly predicted by the Taylor equation 2, in which case a decreasing flow stress would imply a reduction of dislocation density and thus that and cracks could act as dislocation sinks. This possibility is compatible with the idea that tensile cracks correspond to free surfaces within the material, and dislocations intersecting those free surfaces would form steps and disappear from the crystals. One other possibility is that cracks relax internal stress and strain incompatibilities between grains, i.e., act as "geometrically necessary" structures.
A third option is that deformation at low pressure, at which cracks are pervasive, is not controlled by dislocation motion but dominated by elastically-accommodated intergranular slip, and tensile cracks relax the associated internal stresses.
The origin of microcracks during semi-brittle deformation of calcite is potentially coupled to intracrystalline plasticity. As discussed extensively by Nicolas et al. (2017), cracks can be nucleated due to stress concentrations at the head of dislocation pile-ups (e.g., Stroh, 1954;Olsson & Peng, 1976;Wong, 1990) and may therefore be dislocation sinks, i.e., microcracks could contribute to dislocation escape from the deformed crystals. Where slip transfer is inefficient (Olsson & Peng, 1976) and also where geometry results in a high density of GNDs (e.g., Figure 8d) stresses are high, which can lead to the nucleation of cracks.
An approach based on the Kocks-Mecking model has the advantage that twinning and grain size effects can be incorporated, whilst introducing a confining-pressure dependence due to the propagation of Mode I cracks. At this stage however, there are insufficient data on dislocation-density evolution during semi-brittle deformation to make sensible progress beyond the qualitative statements listed above. Thus, detailed modelling attempts are beyond the scope of the present work.

Conclusions
We performed a series of triaxial deformation experiments using three calcite rocks of variable grain size: Solnhofen limestone, Carrara marble and Wombeyan marble. Our experimental results demonstrate that the strength and hardening rates of calcite rocks deformed in the semi-brittle regime are inversely dependent on grain size.
Microstructural observations using visible-light and electron microscopy demonstrate that strain is accommodated by cracking, twinning and dislocation glide. Deformation tests were accompanied by in-situ measurements of P-wave speed, which generally decreases with strain. Wave speed decreases more at room temperature than at 200 • C and 400 • C, which is consistent with microcracking being more prevalent at low temperature.
Quantitative microstructural observations reveal that microcrack density is inversely proportional to hardening rate. While the exact role of cracks in the overall stress-strain behaviour remains unclear at this stage, we propose the hypothesis that tensile microcracks cause weakening by either relaxing internal stresses (accommodating strain incompatilibities and reducing the need for geometrically necessary dislocations), or offering free surfaces where dislocations can escape individual grains, or a combination thereof.
Furthermore, significant decreases in wave speed are observed upon removal of confining pressure, indicating the accumulation of sample damage. Decompression induced wave-speed decreases are greatest in experiments performed at high pressure (600 MPa) and low temperature (room temperature) in conjunction with the highest hardening rates.
Microstructural observations from samples deformed at these conditions are characterised by open grain boundaries, suggesting that wave-speed decreases during decompression originate from the release of stored elastic strain.
-46-manuscript submitted to JGR: Solid Earth Electron microscopy shows that twin density is high, consistent with previous studies at similar conditions. Twin spacing is always at least one order of magnitude smaller than grain size. The spatial distributions of intragranular misorientation suggest that twin boundaries do not always act as significant barriers to dislocation motion and sliptransfer computations indicate that twins are statistically weaker barriers than grain boundaries. Therefore, grain size exerts a first-order control on strength and strain hardening, whereas the spacing of twin boundaries may exert a second-order control on these properties.
Taken together, our results show that semi-brittle flow in calcite is controlled by grain-size dependent processes that lead to significant hardening. This behaviour is qualitatively consistent with rheological models that include dislocation density and twin spacing as key state variables. The role of cracking in the decrease of strain hardening at low pressure and temperature requires the addition of a quantity describing crack density (that should include information on crack spacing, length, and orientation distribution) as a new state variable, to address fully the stress-strain behaviour of rocks in the semibrittle regime.

Acknowledgments
Extensive discussions with Erik Rybacki and Brian Evans, who shared some of their (