SUMMARY
Under compressive stress, rock ‘damage’ in the form of tensile microcracks is coupled to internal slip on microscopic interfaces, such as pre-existing cracks and grain boundaries. In order to characterize the contribution of slip to the overall damage process, we conduct triaxial cyclic loading experiments on Westerly granite, and monitor volumetric strain and elastic wave velocity and anisotropy. Cyclic loading tests show large hysteresis in axial stress–strain behaviour that can be explained entirely by slip. Elastic wave velocity variations are observed only past a yield point, and show hysteresis with incomplete reversibility upon unloading. Irrecoverable volumetric strain and elastic wave velocity drop and anisotropy increase with increasing maximum stress, are amplified during hydrostatic decompression, and decrease logarithmically with time during hydrostatic hold periods after deformation cycles. The mechanical data and change in elastic properties are used to determine the proportion of mechanical work required to generate tensile cracks, which increases as the rock approaches failure but remains small, up to around 10 per cent of the net dissipated work per cycle. The pre-rupture deformation behaviour of rocks is qualitatively compatible with the mechanics of wing cracks. While tensile cracks are the source of large changes in rock physical properties, they are not systematically associated with significant energy dissipation and their aperture and growth is primarily controlled by friction, which exerts a dominant control on rock rheology in the brittle regime. Time-dependent friction along pre-existing shear interfaces explains how tensile cracks can close under static conditions and produce recovery of elastic wave velocities over time.
1 INTRODUCTION
Brittle deformation of rocks is characterized by macroscopic failure and fault slip, and is driven by frictional slip along pre-existing interfaces and tensile fracturing at the microscale. Tensile microfracturing is what is generically termed ‘damage’, and is a key geological marker of brittle deformation in the field. It has been extensively documented in regions surrounding major crustal faults, termed ‘fault damage zones’, both in the field (e.g. Faulkner et al. 2010, references therein) and in situ by seismological methods (e.g. Li et al. 1990, 2004; Ben-Zion & Malin 1991; Allam et al. 2014; Wang et al. 2019). Understanding how ‘damage’ evolves in crustal rocks is crucial for our understanding of the seismic cycle, because the presence of tensile microcracks has a first order impact on both the post-yield rheology and the effective properties of rocks: it increases their porosity, permeability and their elastic compliance and leads to anisotropy when the crack population has a preferred orientation (e.g. Paterson & Wong 2005, chap. 5). In triaxial conditions, stress-induced tensile microcracks tend to be oriented parallel to the most compressive axis, which leads to vertically transverse isotropy in the crack orientation distribution and thus in elastic moduli.
Damage is often characterized and quantified by a measure of density and orientation of open microcracks: this is how fault damage zones have been documented in field and experimental studies (e.g. Chester & Logan 1986; Moore & Lockner 1995; Vermilye & Scholz 1998; Zang et al. 2000; Faulkner et al. 2006; Mitchell & Faulkner 2009; Aben et al. 2020b). Under compressive non-hydrostatic stress states, open microcracks result from local tensile stresses between grains of different elastic properties, at the surface of open pores or at the edges of sliding defects such as grain boundaries (Tapponnier & Brace 1976; Kranz 1983). They are a consequence of shear deformation, but they contribute mostly to volumetric strain, by producing dilation (Brace & Byerlee 1966), and to an increase in compliance perpendicular to their orientation. These dilational effects can be detected seismically (e.g. Gupta 1973; Bonner 1974; Lockner et al. 1977; Schubnel et al. 2003; Aben et al. 2019). By contrast, microscale slip along pre-existing flaws or grain boundaries produces additional shear compliance with hysteresis (e.g. Kachanov 1982a), as documented under uniaxial (David et al. 2012) and triaxial conditions (David et al. 2020), but does not leave any direct microstructural evidence or seismic signature, and is therefore not usually considered as ‘damage’ despite its major role in the micromechanics of brittle rock deformation.
The coupling between internal slip and tensile microcracking is the fundamental basis of our understanding of the rheology and physical properties of crystalline, non-porous rocks in the brittle regime under compressive loading. Microstructural evidence for such coupling was established by Tapponnier & Brace (1976), and many micromechanical models have successfully reproduced the key features of brittle deformation, including stress–strain behaviour, dilatancy, permeability and wave velocity evolution, based on the so-called ‘wing crack’ model, whereby slip on shear cracks is coupled to tensile opening of ‘wings’ at their tips (e.g. Nemat-Nasser & Horii 1982; Horii & Nemat-Nasser 1985; Nemat-Nasser & Obata 1988; Ashby & Sammis 1990; Basista & Gross 1998; Simpson et al. 2001; Deshpande & Evans 2008; Bhat et al. 2011). Spatial and temporal co-location of shear and tensile internal strain has been directly evidenced in laboratory tests by Renard et al. (2019), giving support to the strong coupling arising from fracture mechanics models. Therefore, ‘damage’ as measured in the field is expected to be associated to significant internal slip if it results from applied shear stresses. By contrast, tensile microfracturing can be generated pervasively under dynamic tensile loading, for instance near a fault during earthquake propagation (Doan & Gary 2009). Thus, the partitioning of inelastic deformation between open cracks versus shear cracks is a potential signature of the process at the origin of damage.
The state of damage at depth in the crust is however difficult to assess, due to (1) the strong load-path dependency and (2) the uncertainty regarding long-term persistence of both internal slip and tensile fracturing. External and internal stresses lead to closure or opening of existing (or newly formed) tensile microcracks, so that damage-induced property variations are partially reversible. A significant source of complexity is due to the microscale frictional rheology of shear cracks, which leads to hysteresis in stress–strain behaviour (Zoback & Byerlee 1975b) as well as in dilatancy (Scholz & Kranz 1974; Hadley 1976), wave velocity (Holcomb 1981; Passelègue et al. 2018) and permeability (Zoback & Byerlee 1975a; Mitchell & Faulkner 2008). In addition, the time-dependency of friction leads to delayed crack closure, which results in time-dependent dilatancy and wave velocity recovery (Scholz & Kranz 1974; Holcomb 1981; Brantut 2015; Meyer et al. 2021). It is likely that such time-dependent effects are responsible for at least part of the observed post-seismic variations in wave velocities around faults, which have been documented in many geological settings (e.g. Li et al. 2003; Schaff & Beroza 2004).
In order to understand the impact and possible seismic signature of damage at depth in the crust, the respective contributions of tensile and shear cracks to the overall inelastic behaviour, as well as the reversibility and long-term evolution of the microstructural state, need to be quantified. Here, we follow the approach of Holcomb (1981) and conduct cyclic loading tests in granite under triaxial compression, combining stress–strain measurements with wave velocity monitoring. A similar experimental setup was used recently by Passelègue et al. (2018), who documented in detail the reversibility of stress-induced anisotropy in granite. To go beyond the qualitative description of memory effects as originally reported by Holcomb (1981), we use wave velocities to compute effective moduli of the rock, which are sensitive to open cracks. By subtracting the contribution of those open cracks in the mechanical behaviour, we access separately to the contribution of internal slip. We estimate the partitioning of inelastic strain between microslip and opening as a function of proximity to failure (with increasing stress induced damage). In addition, we quantify the reversibility of dilatancy and elastic wave velocities upon unloading and as a function of time.
2 EXPERIMENTAL METHODS
2.1 Mechanical testing
Samples of Westerly granite were cyclically deformed in a triaxial apparatus (Eccles et al. 2005) at confining pressures of 40, 80 and 120 MPa and with increasing differential stress, under dry conditions. The confining medium was silicone oil, and the confining pressure was controlled within ±0.5 MPa. Axial load was measured externally by a load cell, and corrected for seal friction as in David et al. (2020). Sample shortening was measured with external Linear Variable Differential Transformers (LVDTs), corrected from elastic distortion of the loading column. Stresses and strains are taken positive in compression. The samples were instrumented with two pairs of axial and radial strain gauges, and jacketed in a nitrile sleeve fitted with an array of 14 P- and 2 Sh-sensitive piezoelectric transducers in direct contact with the rock surface. At regular time intervals throughout the tests, the transducer array was used to measure time of flight of P and Sh elastic waves at four different orientations with respect to the compression axis. Only direct arrivals are considered, and converted phases are not analysed. Details of the experimental technique and data processing can be found in Brantut (2015).
The samples were first deformed in load-unload cycles with maximum differential stress gradually increasing by around 100 MPa between successive cycles. In each cycle, a constant axial deformation rate of 10−5 s−1 was imposed until the target differential stress was reached. The samples were then unloaded at the same deformation rate down to minimum differential stress of around 20 MPa. These deformation cycles were conducted at increasing differential stress up to a maximum that was lower than the peak stress (failure strength) of the rock. Then, a second sequence of loading cycles was imposed, with maximum differential stress gradually decreasing by around 200 MPa between successive cycles.
An additional experiment was performed in a similar fashion, except that the load was fully removed after each unloading step, and the sample was held under constant hydrostatic pressure during a 24-hr period before the next loading step.
2.2 Dynamic moduli inversion
During triaxial compression, the elastic properties of rocks become strongly anisotropic, typically developing vertically transverse isotropy (e.g. Schubnel et al. 2003; Passelègue et al. 2018). In anisotropic materials, the times of flight obtained by pulse transmission between transducer pairs located at an angle from the plane of isotropy do not purely correspond to pure P or S phases. Assuming homogeneity, the wave velocities computed from such times of flight are group velocities (e.g. Thomsen 1986). The computation of the stiffness matrix from group wave velocities is a non-linear inverse problem where unknowns also include the phase angles between transducer pairs (Kovalyshen et al. 2017). We solve this problem in the least-square sense using the quasi-Newton algorithm Tarantola (2005, section 3.4.4), taking advantage of automatic differentiation for the computation of the Jacobian matrix required at each iteration.
3 HYSTERESIS IN STRAIN AND ELASTIC WAVE VELOCITIES
Under purely hydrostatic conditions, Westerly granite exhibits some hysteretical behaviour in terms of elastic wave velocities and volumetric strain (Fig. 1): upon a full pressurization–depressurization cycle between 10 and 120 MPa, a permanent increase in VP and VS of around 3 per cent is observed, together with a small permanent decrease in bulk volume, of the order of 0.02 per cent. The hysteresis in wave velocity and strain indicates that some irreversible microcrack closure occurs during hydrostatic loading (Walsh 1965).
Average P- and S-wave velocities and volumetric strain during a hydrostatic pressure cycle.
Under triaxial stress conditions, three main regimes can be distinguished depending on the maximum stress achieved during the cycle and the loading history (Fig. 2). In successive loading cycles at relatively low stress conditions (typically below the onset of inelastic dilation), significant hysteresis was observed in differential stress versus axial strain behaviour, whereas volumetric strain was essentially fully reversible. The elastic wave velocities varied only slightly: with increasing differential stress, horizontal P- and S-wave velocities remained approximately constant and a slight decrease initiated at the highest loads. Subvertical P-wave velocity increased at low differential stress. Some hysteresis also occurs in wave velocities, with values remaining lower during unloading than during loading. All the measured quantities (stress, strain and wave velocities) at the maximum differential stress achieved in a given cycle were matched again when stress returned to the same level in the subsequent cycle.
Evolution of differential stress, axial and volumetric strain and P- and S-wave velocities during representative pairs of consecutive load-unload cycles in Westerly granite at 120 MPa confining pressure. The estimated peak stress at this pressure is (eq. 1) 805 MPa. In left-hand panels, red curves mark the stress–strain behaviour during the first of the two cycles. Legends indicate angles of propagation with respect to axial direction; 90° is perpendicular to loading axis (horizontal) and 0° is along loading axis (vertical).
During loading cycles conducted to higher maximum stress, beyond the onset of measurable dilatancy, hysteresis was observed in both axial and volumetric strain. The horizontal wave velocities evolved in three steps: they remained approximately constant up to a stress threshold (around 300 MPa in the test at Pc = 120 MPa, independently from the cycles), then decreased linearly up to the previous maximum stress experienced by the rock, and then decreased significantly more beyond that stress. Wave velocities changed only very slightly at the onset of unloading, and returned gradually to their initial values with further unloading. The wave velocity recovery was almost total after unloading [as documented by Passelègue et al. (2018)], except for a small permanent decrease. Upon reloading after a given cycle, the axial and volumetric strain as well as the wave velocities followed a qualitatively similar change compared to that in the previous cycle, but did not overlap until the previous maximum stress is reached.
When loading cycles were conducted to a lower maximum differential stress compared to overall maximum achieved during the experiment, the axial strain showed hysteresis, and was fully reversible when load was removed (Fig. 2, bottom left-hand panel). The volumetric strain showed some non-linearity with stress but no hysteresis. Similarly to the pre-dilatancy behaviour, horizontal wave velocities tended to decrease with increasing load. As documented by Holcomb (1981), axial and volumetric strain, and wave velocities completely overlapped during reloading as long as the maximum stress did not exceed the overall maximum achieved during the test.
The stress–strain behaviour and elastic wave velocities of Westerly granite measured here are entirely consistent with previous observations by Zoback & Byerlee (1975b), Holcomb (1981) and Passelègue et al. (2018). Turning our attention to the state of the rock after each unloading stage, we observe a gradual increase in irrecoverable axial and volumetric strain with increasing maximum differential stress (Fig. 3). Irrecoverable dilatant volumetric strain was markedly larger at 40 MPa confining pressure, whereas it remained negligible at 120 MPa. By contrast, more axial permanent strain remained after unloading at high confining pressure compared to low confining pressure. In subsequent cycles conducted to lower maximum stress, permanent strain remained approximately constant.
Irrecoverable axial and volumetric strain after unloading as a function of maximum differential stress achieved during each cycle, Qmax, normalized by the estimated failure stress Qpeak (eq. 1). Black symbols: loading cycles conducted at increasing maximum differential stresses. Grey symbols: loading cycles conducted at decreasing maximum differential stress.
Elastic wave velocities decreased dramatically with increasing differential stress and recovered during unloading, but not always completely (Fig. 4). Permanent changes remained after cycles where the maximum differential stress exceeded about 50 per cent of the estimated reference peak stress, which roughly coincides with the stress threshold for dilatancy (even though such a quantity is difficult to quantify precisely, as stated by Hadley (1975a, p. 22). The permanent drop in wave velocity increased with increasing maximum differential stress, and remained constant in any subsequent cycles where the previously achieved maximum differential stress was not overcome. Except for one cycle at high differential stress at 40 MPa confining pressure, residual wave velocity drops were typically of the order of a few percent. The elastic anisotropy that developed at high stress weakened but also remained non-zero after unloading (keeping in mind the caveat that unloading was not conducted down to zero but to about 20 MPa differential stress), with subvertical P-wave velocities being higher than subhorizontal ones. The residual P-wave anisotropy can be characterized by Thomsen’s ϵ parameter defined as ϵ = (C11 − C33)/(2C33) (where Cij denote components of the elastic tensor in Voigt’s notation), which showed a decrease towards more negative values after stress cycles with increasing differential stress (Fig. 5). At 40 MPa confining pressure, a relatively large permanent anisotropy remained (ϵ ≈ −0.125), in conjunction with a dilatant residual volumetric strain. With increasing confining pressure, the residual anisotropy decreased, down to around ϵ = −0.025 after the largest stress cycle at 120 MPa confining pressure. In all cases, the anisotropy did not change significantly with subsequent loading cycles when the maximum stress achieved was less than the overall maximum reached during the test. The residual anisotropy was correlated to the residual volumetric strain (Fig. 5).
Relative change in elastic wave velocities at the maximum differential stress (top panels) and after unloading (bottom panels).
Thomsen’s ϵ parameter measured after unloading as a function of (a) maximum differential stress achieved during each cycle, and (b) irrecoverable volumetric strain. Black symbols: loading cycles conducted at increasing maximum differential stresses. Grey symbols: loading cycles conducted at decreasing maximum differential stress.
After each experiment, elastic wave velocities decreased during hydrostatic decompression, down to values lower than those of the intact material (e.g. by about 10 per cent for the test conducted at 120 MPa confining pressure, cf. Fig. 6). The decrease in elastic wave velocities was accompanied by a significant increase in anisotropy (Fig. 7).
Change in elastic wave velocities during hydrostatic decompression after a series of 16 loading cycles up to 840 MPa differential stress.
Evolution of Thomsen’s parameter ϵ during hydrostatic decompression after deformation.
4 ROLE OF TENSILE MICROCRACKS
Thin sections were made from an intact sample and from the sample deformed at 80 MPa confining pressure, which experienced a maximum differential stress of 721 MPa and a total of 14 load-unload cycles. Despite a post-deformation wave velocity at least 20 per cent lower than the intact rock, qualitative optical observations reveal very little, if any, signs of mode I intragranular microcrack damage beyond what is observed in the intact material (Fig. 8). While it is well known that detection of stress-induced intragranular cracks is difficult with basic optical methods (Tapponnier & Brace 1976), the apparent discrepancy between the large change in properties and the high stress level experienced by the sample and the lack of any widespread, visible microstructural evidence is quite striking. That being said, our observations are still consistent with the detailed microstructural work of Tapponnier & Brace (1976), who could only detect newly formed intragranular cracks at stresses significantly higher than the onset of dilatancy.
Optical micrographs of (a) intact and (b) deformed Westerly granite, in transmitted plane polarized light. Compression axis is vertical.
Therefore, it is likely that a significant source of dilatancy and wave velocity variations observed during the test, especially at low stress, was the opening of pre-existing cracks or grain boundaries. The large recovery of properties during unloading implies that any newly open crack could be almost perfectly closable, that is crack faces should return into contact (at least partially) as load is removed.
Assuming that the anisotropic variations in elastic moduli arise from open microcracks, we use an effective medium model to invert the time of flight data for quantitative crack density parameters. Here we follow the approach of Sayers & Kachanov (1995), who consider an array of uniformly distributed penny-shaped microcracks embedded in a homogeneous, isotropic elastic matrix. Interactions between cracks are neglected, which is valid in the limit of small crack density, and crack orientations are assumed to be distributed in a transversely isotropic manner, with symmetry axis along the most compressive load direction. Sayers & Kachanov (1995) express the effective elastic properties of the cracked material in terms of the moduli of solid and five independent crack density parameters: two components α11 and α33 of a second-order crack density tensor, and three components β1111, β1133 and β3333 of a fourth-order crack density tensor. The components α11 and α33 correspond to the vertical and horizontal crack densities defined by Nvc3, where Nv is the number of cracks per unit volume and c is their average radius. The components βijkl correspond to fourth order moments of the crack orientation distribution. Complete mathematical definitions can be found in Sayers & Kachanov (1995). The solid matrix Young’s modulus is taken equal to E = 89 GPa and its Poisson’s ratio is ν = 0.22, chosen to match intact P- and S-wave velocities of 6.2 and 3.7 km s–1, respectively. To better constrain the inversion, we use the prior knowledge that the β’s are a factor ν/2 smaller than the α’s (Sayers & Kachanov 1995, their eq. 15), that is roughly one order of magnitude smaller. The inverse problem is solved by the quasi-Newton method of Tarantola (2005, section 3.4.4), using the time of flights as observables and crack density parameters and phase angles as unknowns. Since we expect the α’s to be of the order of 0.1, we use a prior value of 0 and a variance of 10−2 on all β’s.
The evolution of horizontal and vertical crack densities mirrors the evolution in elastic wave velocities (illustrated for representative cycles at 120 MPa confining pressure in Fig. 9). In the low stress regime, the horizontal crack density decreases slightly with increasing stress, while the vertical crack density remains stable. Beyond the onset of dilatancy, the vertical crack density starts increasing and the horizontal crack density tends to stabilize. For cycles conducted up to relatively low differential stress, hysteresis is observed only in the vertical crack density, which is almost completely recoverable upon full unloading. At elevated differential stress, the vertical crack density increases with increasing stress, until the previous maximum stress is reached, and increases more markedly with further loading. Upon unloading, the vertical crack density remains high, and recovers only partially. The horizontal crack density remains more or less stable at high stress, and increases linearly with differential stress during unloading, showing hysteresis. In cycles conducted to lower maximum differential stresses compared to previous cycles, the evolution of crack density is perfectly reproducible during reloading. No hysteresis is observed in horizontal crack density, but hysteresis remains in the evolution of vertical crack density.
Evolution of vertical (α11) and horizontal (α33) crack density with differential stress during representative pairs of consecutive load-unload cycles in Westerly granite at 120 MPa confining pressure. Red curves mark the behaviour during the first of the two cycles.
Overall, the density of open, mostly vertical, microcracks evolve non-linearly with stress and decrease (i.e. open microcracks mechanically close) with strong hysteresis during unloading. Pre-existing cracks make the material behave non-linearly, that is compliance decreases with increasing stress. The formation of new cracks, either by extension or re-opening of pre-existing ones, or by creation of new crack surfaces in intact material, is achieved only when the previously achieved maximum stress is overcome. When no new cracks are generated in a given load cycle, the crack density evolution during loading is reproducible, but hysteresis during unloading remains.
The residual crack density after unloading increases with increasing differential stress, in proportion to the maximum crack density reached during each cycle (Fig. 10). The maximum vertical crack density is of the order of 0.1 when differential stress approaches the rock strength at each confining pressure investigated; this value is consistent with the non-interaction approximation of Sayers & Kachanov’s effective medium model. The residual vertical crack density can be as high as 0.025 after unloading, and remains constant with any further cycles if the previously achieved maximum stress is not overcome.
Vertical (α11) horizontal (α33) crack densities at the maximum differential stress (top panels) and after unloading (bottom panels) for each loading cycle. Black symbols: loading cycles conducted at increasing maximum differential stresses. Grey symbols: loading cycles conducted at decreasing maximum differential stress.
During hydrostatic decompression, both vertical and horizontal crack densities increase with decreasing confining pressure, reaching values of the order of 0.1 at the lowest pressure (Fig. 11). The increase in anisotropy noted in Fig. 7 is explained by a stronger increase in vertical compared to horizontal crack density.
Vertical (black symbols) and horizontal (grey symbols) crack density during hydrostatic decompression after deformation.
5 STRAIN AND ENERGY PARTITIONING: INELASTIC VERSUS ELASTIC EFFECTS OF MICROCRACKS AND MICROSCALE FRICTION
At all the tested confining pressures, the axial strain is dominated by the elastic contribution together with inelastic processes (Fig. 12). The contribution of crack-induced compliance changes is practically negligible. In terms of volumetric strain, the elastic contribution is significant, producing additional compression and the contribution of inelastic processes increases with increasing differential stress to produce dilation, which eventually dominates the total response. There again, purely elastic variations in compliance due to open cracks do not produce significant strains. In any case, the effect of open cracks is to increase elastic compliances in dry rocks, which would tend to produce more compressive strains with increasing damage. The large dilatancy is therefore an effect of inelastic variations in crack length, sliding and crack opening.
Axial (top panels) and volumetric (bottom panels) strain at the maximum differential stress of each cycle. Dashed lines are the elastic contributions of the solid, red curves at the elastic contributions of open microcracks, and black curves are the remaining inelastic contributions, which include internal slip and irreversible crack growth. All contributions sum to produce the observed stress–strain curves. The elastic contribution of open cracks is derived from inversion of wave velocities, and the contribution of slip and crack growth is derived by difference between measured strains and inverted elastic strains from other contributions.
Throughout the loading cycles, the compliance keeps changing, decreasing during unloading and increasing during re-loading (Fig. 2). However, not all these changes correspond to the formation of new cracks: the decrease in compliance during unloading is due to the formation of new contact points along pre-existing cracks (no healing occurs), and the increase in compliance up to the previous maximum load is due to the re-opening of pre-existing cracks. Therefore, only the compliance increase during the loading steps beyond any previously achieved maximum differential stress should be considered in the estimation of fracture energy by eq. (4) (see path BC in Fig. A1). Despite this precaution, it is quite clear that more compliance variations would also be due to re-opening (i.e. removal of contact points) of pre-existing cracks at all stages of the loading process, so that our estimate of fracture energy is an upper bound. An alternative method to obtain a less severe upper bound on integrated fracture energy per cycle from stress–strain curves is discussed in the Appendix.
(a) Energy dissipation per cycle. Black curves correspond to the total dissipation, and red curves correspond to the contribution of microcrack growth computed from eq. (4), which is equivalent to the total fracture energy per unit volume. (b) Ratio of cumulated fracture energy over total dissipation per unit volume as a function of differential stress during loading stages only.
A significant part of the total dissipation by frictional processes occurs during the unloading and reloading steps of each cycle, whereas crack growth occurs when load is increased beyond the previously achieved maximum (segment BC in Fig. A1). The energy budget computed during only those loading steps indicates a growing contribution of fracturing with increasing stress, from around 10 to 40 per cent of the total dissipation near the peak stress (Fig. 13b). At elevated confining pressure, the cumulated fracture energy (per unit rock volume) is a larger fraction of the total dissipation, whereas the inferred crack densities are comparable or even lower than at low pressure (Fig. 10). This indicates that deformation promotes the growth of fewer, larger cracks at low pressure and more numerous, shorter cracks at elevated pressure.
6 TIME DEPENDENCY AND RECOVERY
In order to investigate the long-term persistence of the residual strains and velocity drops that remain after deformation and unloading, one experiment was conducted at Pc = 80 MPa where the rock was held under hydrostatic conditions (with axial load completely removed) for 24 hr between each loading cycle. Elastic wave velocities tend to recover during the hold period, by only about 0.1–0.2 per cent after load cycles up to modest differential stress (Qmax/Qpeak < 0.8), and by up to 0.8 per cent when the maximum stress was near the failure strength (Fig. 14). The recovery is logarithmic with time, and is more marked along subhorizontal orientations. The recovery in wave velocity is small (typically a fraction of a percent) in relation to the wave velocity of the intact material, but represents a large fraction of the net drop that occurred during each cycle: for instance, in the last load cycle, the horizontal P-wave velocity dropped by around 110 m s–1 (residual after unloading), and recovered by around 50 m s–1 over 24 hr, that is a 45 per cent recovery rate.
Relative change in elastic wave velocities as a function of time during hydrostatic hold periods after each loading cycles. Test conducted at 80 MPa confining pressure.
The recovery in wave velocities can be interpreted as a logarithmic decrease in crack density, mostly in the axial direction (Fig. 15a). It is accompanied with a small but detectable logarithmic compactive change in volumetric strain, together with a decrease in elastic anisotropy (Figs 15b and c). The recovery is almost negligible when the maximum differential stress was less than around 70 per cent of the strength of the rock: this is expected since only very small residual change in wave velocity and volumetric strain occur after those load cycles.
Evolution of vertical crack density, volumetric strain and Thomsen’s ϵ parameter as a function of time during hydrostatic hold periods after each loading cycle. Test conducted at 80 MPa confining pressure.
The recovery data are consistent with the time-dependent dilatancy recovery observed by Scholz & Kranz (1974) and Holcomb (1981) in dry granite, and with anelastic strain recovery observed by Gao et al. (2014) in granite, marble and sandstone. In addition, our elastic wave velocity data follow the same qualitative trends as those observed in porous limestone (Brantut 2015) and dry Carrara marble (Meyer et al. 2021).
7 DISCUSSION
7.1 Summary of results
The main experimental observations can be summarized as follows. During deformation at relatively low stress, elastic wave velocities change only slightly, volumetric strain is reversible without hysteresis and shear strain is largely reversible but significant hysteresis is observed. When deformation is continued up to elevated stress, typically beyond the onset of dilatancy, both volumetric and shear deformation, as well as elastic wave velocities, exhibit hysteresis. When load is removed, small residual permanent changes exist (net dilation, shear strain and velocity drop), and recover logarithmically with time under hydrostatic conditions. When the samples are re-loaded up to differential stresses lower than any previous maximum, shear strain and elastic wave velocities exhibit hysteresis, but no further permanent change is recorded and the reloading paths all overlap. The small irrecoverable changes in elastic wave velocity (and associated anisotropy) become larger during hydrostatic decompression.
The stress–strain behaviour observed here is identical to the results of Zoback & Byerlee (1975b), and their interpretation in terms of coupled grain boundary friction and microcrack opening (a qualitative description of what became later the ‘wing crack’ model, originated by Brace et al. (1966)) remains largely valid in the light of complementary measurements of elastic wave velocities, including those from previous authors (Holcomb 1981; Passelègue et al. 2018). During loading, the decrease in elastic moduli is interpreted as the opening microcracks preferentially oriented along the compression axis (consistent with the anisotropic decrease in wave velocities), while inelastic axial (shear) strain may be attributed to slip along pre-existing interfaces, most likely grain boundaries.
Slip along pre-existing flaws is not the only mechanism that can generate tensile microcracks: stress concentrations due to open pores (e.g. Sammis & Ashby 1986) and elastic anisotropy or moduli mismatch at grain contacts (e.g. Dey & Wang 1981) also lead to crack nucleation and growth. In Westerly granite, there is no detectable equant porosity, so pore-emanating cracks are not significant. Tensile cracking due to moduli mismatch was analysed in a simple configuration by Dey & Wang (1981); in this model, tensile cracking is linked to shear stresses that are proportional to differences in elastic moduli between contacting grains, but is not linked to actual slip at the contact. From a kinematic point of view, axial tensile cracks originating from this mechanism are not associated to additional axial strain and are not expected to produce hysteresis. The model of Dey & Wang (1981) relies on a number of simplifying assumptions, and further analysis would be needed to determine how shear and tensile cracks are coupled when elastic mismatch is significant. Nevertheless, while this mechanism is potentially relevant to explain tensile cracking in our experiment on Westerly granite, it is insufficient to explain the observed hysteresis in stress–strain curves.
While the role of internal slip in the generation of tensile microcracks has been well established for several decades (see review in Paterson & Wong 2005, section 5.8), a number of features were still undocumented. In particular, the energy budget analysis shows that frictional processes largely dominate the stress–strain response, and the contribution of crack propagation becomes non-negligible (but remains quite small) only at the highest loads. In addition, our results confirm those of Holcomb (1981); Passelègue et al. (2018) in that most of the elastic wave velocity variations are recoverable upon unloading, but also indicate the existence of small but measurable irrecoverable changes in elastic properties, volumetric strain and anisotropy. The residual changes increase with increasing maximum stress, and are amplified during hydrostatic decompression.
7.2 Role of friction
In order to analyse the role of friction on the deformation and physical properties of brittle rocks, it is helpful to use a micromechanical framework that contains the key necessary ingredients: slip on pre-existing grain boundaries or defects and opening of tensile cracks. A number of micromechanical models of brittle deformation have been published that investigate the coupled evolution of slip and crack opening and its effect on macroscopic strain. Perhaps the first model that accounted for friction along pre-existing cracks is that of McClintock & Walsh (1962), who computed growth criteria for such cracks; the effect of friction has subsequently been included in many models aimed to determine thresholds for inelastic behaviour or dilatancy (Nemat-Nasser & Horii 1982; Horii & Nemat-Nasser 1985; Ashby & Hallam 1986; Ashby & Sammis 1990, among others). The effect of friction on stress–strain behaviour was first analysed quantitatively by Holcomb (1978), with a focus on explaining dilatancy. Further work systematically investigated the coupled role of frictional slip and tensile crack growth on full stress–strain curves, using a variety of approximations to solve the crack problems and volume averaging to obtain macroscopic stress–strain curves (or tangent moduli). Kachanov (1982a, b) solved 3-D crack sliding and extension problems, and used non-interaction approximation to average the strain contributions of cracks to total strain. A similar approach with a simplified crack geometry was used by Fanella & Krajcinovic (1988). Horii & Nemat-Nasser (1983) used a self-consistent effective medium approach to approximate crack interactions and obtain tangent moduli of rocks containing open and sliding cracks. Complete stress–strain curves were obtained by (Moss & Gupta 1982; Nemat-Nasser & Obata 1988) based on a 2-D sliding crack model including dilation. Similar results were obtained by Basista & Gross (1998), who used the energy approach of Rice (1975) to compute macroscopic strains from microscale processes. Similar energy-based micromechanical models were extended to 3-D cases by Pensée et al. (2002), Deshpande & Evans (2008) and Bhat et al. (2011). Overall, all these models vary in their quantitative predictions due to a number of detailed assumptions about crack geometry, kinetic relations and homogenization procedures, but share the same qualitative features.
The physical intuition provided by eq. (6) in a simple model is entirely consistent with our observation of stress–strain behaviour, which show that internal slip dominates the macroscopic response (Figs 12 and13). The amount of internal slip needed to explain the inelastic strain observed during the experiments can be roughly quantified as |$b/c \approx \epsilon ^\mathrm{i}_{11}/(\omega _0 \sin 2\phi )$|. Using a crack density of the order of 0.2, ϕ = π/4 and a maximum inelastic axial strain of the order of 0.5 per cent, the ratio b/c is of around 0.025. Using an initial shear crack size of the order of the grain size, 100 μm, an upper bound for the finite slip at the maximum differential stress is of 2.5 μm. Upon unloading, the residual, irrecoverable axial strain is of the order of 0.15 per cent (Fig. 3), so that the residual internal slip is of less than 1 μm. Such small amounts of slip are unlikely to be detectable in the microstructure, even more so than microcracks re-open during decompression and might mask small offsets between grains or along pre-existing flaws.
The coupling between tensile crack size and slip is apparent in the second term of eq. (6): tensile ‘wings’ provide an additional contribution of slip to lateral strain, that is to dilatancy. This is apparent, for instance, in the correlation between residual volumetric strain and residual anisotropy (Fig. 5): hysteresis due to friction prevents complete closure of tensile cracks, which tend to be mostly in the axial direction, so that dilatant strain evolves in tandem with anisotropy. One can also directly observe that crack growth alone is not expected to produce dilatancy without the coupling with shear stress and slip on the pre-existing flaws. This feature is linked to the rather simple crack geometry of the model with purely axial tensile wings, but again shows the dominant effects.
The micromechanics of deformation arising from the wing-crack model also explains the apparent increase in Poisson’s ratio at high stress, sometimes above 0.5, observed in a number of studies (e.g. Faulkner et al. 2006; Heap et al. 2010). As noted in eq. (6), the ratio of lateral to axial strain increases notably if sliding and crack growth/closure are activated. There is however a very large difference between the tangent moduli (and the apparent Poisson’s ratio that could be computed from stress–strain curves) and the elastic moduli of the cracked rock (elastic Poisson’s ratio of damaged dry rocks typically decrease with increasing crack density). Large Poisson’s ratio can only be explained by inelastic effects (slip and crack growth/closure), which are activated by stress. In other words, large Poisson’s ratio are not an inherent property of ‘damaged’ rocks, but appear at high stresses due to the activation of internal slip. The effect of ‘damage’, or in a narrower sense, tensile microcracks, is to lower the stress threshold for non-linear behaviour, since the fracture toughness that gives rise to cohesion becomes zero when tensile cracks pre-exist. This behaviour is consistent with the ‘erosion’ of the onset of dilation C’ with load cycles under triaxial conditions (Hadley 1976). Therefore, the apparent Poisson’s ratio of rock depends on stress and stress history, and cannot be considered a material constant that depends solely on pre-existing microstructural features.
7.3 Characteristics of microscale friction
During unloading, the inferred open crack density almost immediately decreases as soon as the load is reduced (Fig. 9). As noted by Stevens & Holcomb (1980) and Holcomb (1981), there is apparently no ‘dead band’ in a strict sense, contrary to what wing crack models would predict. However, the initial rate of crack density reduction is initially very small, and increases significantly as load is further reduced. The contrast in behaviour between successive cycles at increasing loads (Fig. 2, cycles 10 and 11), where wave velocities show almost a ‘true’ deadband, and those at decreasing loads (Fig. 2, cycles 13 and 14) where the deadband is less apparent, is compatible with the wing crack concept. During cycles performed up to stresses lower than the previous maximum, the rock has kept the memory of the previously achieved maximum in the form of the elastic restoring force between the shear crack faces, so that reverse sliding can be activated more easily as soon as load is reduced, compared to cases where a new maximum has just been reached. Experimental observations, including those of Holcomb (1981), do not completely invalidate wing crack models, but rather call for a few amendments to the description of wing crack mechanics. Several non-exclusive possibilities exist to explain the near-instantaneous change in crack density upon unloading: (1) complex friction law, for example including direct stress dependency on friction coefficient or hysteresis, (2) viscous relaxation in the bulk (or along shear cracks or grain boundaries) and (3) elastic changes in crack aperture leading to the formation of contact points along crack faces, which have instantaneous stiffening effects on effective elastic moduli (Trofimov et al. 2017; Passelègue et al. 2018; Meyer et al. 2021).
Recovery is also compatible with bulk relaxation and the generation of contact islands along cracks, as discussed in Meyer et al. (2021). In the case of dry granite at room temperature, it is unlikely that viscous flow or fluid-related processes (such as pressure solution) can have any significant role in the recovery phenomenon, and backsliding (accompanied by crack closure) is most likely the dominant process.
The mechanism of backsliding requires that the frictional strength along pre-existing defects such as grain boundaries does not significantly increase as the sense of shear is reversed. In our tests, under dry conditions, at relatively slow rates and with limited net slip (with no direct evidence left in the microstructure), there is no reason to believe that frictional strength can be dramatically different upon reversal of slip direction. This may not be the case in general, for example if forward slip is changing permanently the structure of the interface.
What sort of friction law would be required to explain all experimental observations? First, we should recall the key observation that the dilatancy threshold increases with increasing confining pressure (Brace & Byerlee 1966), which is consistent with the idea that a well-defined friction coefficient (a constant ratio between shear and normal stress at the onset of slip) exists, to first order. The near instantaneous, small backslip upon unloading inferred from the variation of wave velocities is compatible with the existence of an elastic shear stiffness of the microscale interfaces. The time-dependent recovery is compatible with rate-strengthening friction, as described (for instance) by rate-and-state laws. Therefore, friction at the microscale (grain contacts being of the order of 100 μm) has all the characteristics of macroscopic friction. Models of macroscopic friction rely on rough interfaces with elastic or plastic microcontacts (e.g. Scholz 2019, chap. 2), and the microscopic frictional interfaces are therefore expected to be also rough. This is potentially due to the opening of grain boundaries due to cooling and exhumation, forming geometric mismatch between internal interfaces.
We expect microscale friction to change character at high pressure or high temperature, if all open cavities are closed: this would correspond to saturation of contact area, at which point resistance to slip is not directly proportional to normal stress but governed by interfacial defects (e.g. grain boundary dislocations) or adhesive forces. Such a transition has been inferred in experiments conducted in antigorite at high pressure and temperature (David et al. 2020). The change in characteristics of microscale friction may contribute to the brittle–ductile transition, in addition to commonly recognized processes such as the activation of intracrystalline creep mechanisms.
7.4 Contribution of tensile cracking to faulting
Throughout this paper we have considered the evolution of average rock properties with deformation cycles. However, at the highest differential stresses reached in some of the tests, we expect strain to be at least partially localized. Hadley (1975b) showed that circumferential strain loses circular symmetry around the compression axis at differential stress beyond 90 per cent of the failure stress, indicating an early tendency towards the formation of a shear fault. In situ, time-resolved X-ray tomography observations (Renard et al. 2019) seem to indicate that strain localization in the form of shear and tensile crack coalescence is a rather late process, occurring at stresses around 99 per cent of the peak. The inception of macroscopic shear faulting is caused by microcrack interactions, which facilitate further tensile crack propagation and slip on shear cracks. Despite these additional quantitative changes to the mechanics of microcracking and internal slip, interactions do not change the essential qualitative characteristics of friction and tensile crack opening at the microscale.
Damage accumulation (per inelastic strain increment) in the form of tensile microcracks is known to accelerate as macroscopic failure approaches and strain becomes localized (e.g. Scholz 2002, section 1.2.2). At this stage, the contribution of tensile cracking to the overall energy dissipation may deviate from what is observed at lower stress levels. Our data show a trend towards a more prominent contribution of open cracks at high stress (Fig. 13). In addition, the cumulative energy dissipation from microcracks is also substantial, up to 40 per cent or more (if, say, the rock was brought to shear failure in a single cycle). Thus, it is possible that the energy required for tensile crack growth could transiently dominate the energy budget in the phase immediately preceding shear failure, where differential stress is near its peak.
The transition from homegeneous deformation to shear failure cannot be analysed with our laboratory methods, which require volume averaging of strains and elastic wave velocities. Nevertheless, the energy balance of pre-rupture deformation given here (Fig. 13) is consistent with that obtained from P-wave velocity tomography by Aben et al. (2019, 2020a) on faulted rock: during quasi-static shear failure propagation, the contribution of tensile fracturing as captured by changes in elastic moduli remains a small fraction (less than 10 per cent) of the total energy dissipated by other processes such as friction. More tensile crack damage is generated during dynamic faulting than in quasi-static tests, but the change is not dramatic (Aben et al. 2020a). The localized deformation around a shear fault (fault ‘damage’ zone) is therefore of the same nature as the bulk deformation of the nominally intact rock, and is likely dominated by shear.
Tensile cracks do not need to be present in large densities or have pervasively large dimensions throughout the rock volume to produce large effects on crack coalescence, shear failure and slip. This is consistent with the fact that crack coalescence and catastrophic failure is linked to local maxima in stress intensity factors (Kachanov & Sevostianov 2012), that is only the most deleterious cracks play a role.
8 CONCLUSIONS
A dominant component of stress-induced ‘damage’ in rocks can be attributed to microscale slip on pre-existing interfaces, such as grain boundaries, which is largely reversible upon unloading and does not leave any substantial microstructural record. Tensile cracking is coupled to slip, and produces large changes in dynamic elastic properties and dilatancy, but only contributes to a small fraction of the energy required for deformation. As originally discussed by Holcomb (1981), rocks keep a memory of the maximum differential stress applied to them. The simplest conceptual model to explain this memory effect is that the maximum stress is recorded as an elastic restoring force along frictional interfaces, well captured by wing crack models. Residual, irrecoverable strains, elastic wave velocity drops and anisotropy can be attributed to the hysteresis associated with backslip, where friction has to be overcome in the reverse slip direction. Complete decompression amplifies the small residual changes that occurred during triaxial loading.
The inelastic behaviour governed by friction and tensile fracturing implies that deformation and physical properties depend on stress path and not only on the current stress state. Due to the threshold effect of friction, ‘damage’ makes rock behave in a non-linear fashion, and is only visible at elevated stress (where slip and crack opening can be activated). In the crust, differential stress is limited by macroscopic frictional strength of faults (e.g. Zoback & Healy 1992; Zoback & Harjes 1997), which has the same characteristics as frictional strength of microscale interfaces responsible for inelastic rock behaviour (as long as the fault core composition is similar to that of the wall rock). Therefore, one expects that the manifestation of ‘damage’ can only be significant in regions where local differential stress is (or has been) sufficiently large to allow for internal slip to be activated: ‘damage’ has significant consequences only at elevated differential stress. If a volume of rock has been subjected to transiently high compressive shear stress, it can keep a memory of the peak stress in its elastic anisotropy, but this memory can be easily erased with further changes in stress orientation and magnitude. This has been evidenced in laboratory true triaxial experiments (Browning et al. 2017, 2018). In nature, we expect shear stress to fluctuate cyclically in the vincinity of seismogenic faults. Although shear stress is on average likely close to ‘Byerlee’-type frictional strength, it can be locally much higher (e.g. at geometric irregularities). In addition, during earthquake propagation, near-tip stresses are also transiently well above the long-term applied stress level. Thus, rocks located close to seismogneic faults will experience stress cycles during which microscale slip and tensile cracking will occur. It is in these so-called ‘damage’ zones that significant strain energy is dissipated during the seismic cycle (e.g. Griffith et al. 2010; Okubo et al. 2019), and we highlight here that tensile cracking may not be the dominant contribution to this dissipation.
Time-dependent recovery in elastic wave velocities and dilatancy under constant stress conditions is consistent with rate-dependent friction along microscopic interfaces, coupled to backslip-induced crack closure. Time-dependent friction implies that rocks have a viscous behaviour even under low temperature, low pressure and dry conditions. Detailed modelling and experimental work is required to better characterize this time dependency and how it interacts with known other time-dependent aspects of brittle rock deformation, such as subcritical cracking (Brantut et al. 2013). Nevertheless, the possibility of bulk stress relaxation and inelastic deformation driven by time-dependent microslip has potentially large consequences for our understanding of transient deformation phenomena in the brittle crust such as post-seismic deformation and elastic wave velocity variations during the seismic cycle.
ACKNOWLEDGEMENTS
Frans Aben and Emmanuel David helped running the experiments. Teng-Fong Wong and an anonymous reviewer made useful suggestions that helped clarify the paper. Support from the UK Natural Environment Research Council (grant NE/K009656/1) and from the European Research Council under the European Union’s Horizon 2020 research and innovation programme (project RockDEaF, grant agreement #804685), is gratefully acknowledged.
DATA AVAILABILITY
Data are available at doi:10.5281/zenodo.6761111. The code used to invert group velocity data for elastic moduli is accessible at https://github.com/nbrantut/VTIModuli.jl. The code version used here was commit dae4376.
REFERENCES
APPENDIX: ALTERNATIVE ESTIMATE OF INTEGRATED FRACTURE ENERGY
During a complete loading and unloading cycle, the energy dissipated is the sum of contributions including friction, crack opening and crack growth. Energy changes due to tensile crack growth can be estimated from independent measurements of elastic moduli (Section 5), but are also reflected in the overall stress–strain behaviour. Here, we establish an approximate technique to estimate an upper bound of the energy change due to tensile fracturing during a given loading cycle, based solely on stress–strain analysis.
Let us consider a rock sample that experienced a complete loading cycle (ABCD in Fig. A1), and is now in a macroscopically unloaded state (point D in Fig. A1). Upon reloading from point D, the stress remains lower than during the previous cycle, but meets again point C, the maximum in stress and strain achieved in the previous cycle. The energy dissipated in the unloading-loading cycle CDC is purely frictional, since we do not expect any significant propagation of tensile cracks as long as the applied stress remains lower than the previously achieved maximum. The difference between the total dissipated energy in the cycle (area ABCD) and the frictional dissipation upon reloading (area CDC), denoted Wc, must therefore include all the contribution of tensile crack propagation during that cycle. It certainly includes other contributions too, such as that of forward sliding of cracks during loading from A to B, which occurs at a higher stress state and lower tensile crack density than from C to D. Wc is thus an upper bound for the contribution of tensile crack growth in a given cycle. When computing the fracture energy from changes in elastic moduli ΔG (eq. 4), only the portion BC is considered to minimize the contribution of reversible changes in moduli due to re-opening of pre-existing cracks.
Sketch of stress–strain behaviour during successive loading cycles. The cycle of interest initiates in state A and ending in state D. The red portion BC corresponds to the stress and strain interval where the previous maximum stress is overcome, that is where tensile cracks further propagate and new ones are formed.
Comparison of total energy dissipation per cycle (Wtot), fracture energy estimates ΔG from changes in elastic moduli and Wc from stress–strain analysis.
