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

# Deformation evolves from shear to extensile in rocks due to energy optimization

Jess McBeck<fnref n="1" />, Benoît Cordonnier², Michele Cooke<fnref n="3" />, Laura Fattaruso<fnref n="3" /> & François Renard<fnref n="1" />,⁴

Determining how fracture network development leads to macroscopic failure in heterogeneous materials may help estimate the timing of failure in rocks in the upper crust as well as in engineered structures. The proportion of extensile and shear deformation produced by fracture development indicates the appropriate failure criteria to apply, and thus is a key constraint in such an effort. Here, we measure the volume proportion of extensile and shear fractures using the orientation of the fractures that develop in triaxial compression experiments in which fractures are identified using dynamic in situ synchrotron X-ray imaging. The fracture orientations transition from shear to extensile approaching macroscopic, system-size failure. Numerical models suggest that this transition occurs because the fracture networks evolve in order to optimize the total mechanical efficiency of the system. Our results provide a physical interpretation of the empirical internal friction coefficient in rocks.

<footnote n="1">¹Departments of Geoscience and Physics, Njord Centre, University of Oslo, Oslo, Norway. ²European Synchrotron Radiation Facility, Grenoble, France.</footnote>

<footnote n="3">³Department of Geoscience, University of Massachusetts, Amherst, MA, USA. ⁴ISTerre, Univ. Grenoble Alpes, Grenoble INP, Univ. Savoie Mont Blanc, CNRS, IRD, Univ. Gustave Eiffel, 38000 Grenoble, France. ☎email: j.a.mcbeck@geo.uio.no</footnote>
In the brittle regime, at lower confining stress and depth within the Earth's crust, heterogeneities influence the macroscopic strength of rocks by controlling where and how fractures nucleate, propagate, and subsequently coalesce¹⁻³. Recent work suggests that fracture network development controls the predictability of failure⁴⁻⁶. Moreover, characteristics of fracture networks revealed in X-ray tomography triaxial compression experiments, such as the length and orientation, can predict with success the timing of macroscopic failure and the growth of individual fractures⁷,⁸. Theories from linear elastic fracture mechanics predict with success the conditions for propagation, and geometry of isolated, individual fractures⁹,¹⁰. However, these theories struggle to characterize and predict the coalescence of several thousands of fractures as rocks and other heterogeneous material approach macroscopic failure under compres- sion e.g.,²,¹¹. One of the fundamental characteristics of the deformation produced by fracture development is the failure mode, the endmembers of which are opening (mode I) and shear (mode II/III). Constraining the relative importance of these two modes of deformation during triaxial compression of hetero- geneous materials has been a fundamental goal in fracture mechanics because the evolving competition between these modes determines the applicability of different failure criteria, and the corresponding theories of the micromechanics of brittle deformation in rocks. In particular, whether a numerical model or analytical analysis assumes that fracture propagation occurs due to tensile failure or shear failure determines the predicted geometry of fracture propagation¹¹.

**The competition between opening and shear deformation in rocks.** Some laboratory observations suggest that in triaxial compression with confining stresses lower than about 10 MPa, the fracture network is dominated by extensile fractures that develop parallel to the direction of maximum compression, σ₁, and open perpendicular to this direction throughout loading³. We adopt here the terminology of “extensile fractures” because “extension” refers to strain, whereas “tension” refers to stress e.g.,³. The fracture geometry and resulting fragments of rock following macroscopic failure indicate that axial splitting and thus extensile-dominated deformation is common in uniaxial experiments³. The application of a few megapascals of confining stress suppresses axial splitting, and promotes the development of fracture networks aligned obliquely to the σ₁ direction that accommodate shear³. Laboratory analyses suggest that the system-spanning fracture networks that accommodate shear develop from the coalescence of preexisting extensile fractures¹²⁻¹⁷, and this process may include the buckling and/or rotation of columns between the extensile fractures¹². Others argue that the propagation of a shear fracture produces its own process zone of extensile deformation that includes a diffuse cloud of extensile fractures ahead of the shear fracture¹⁸.

However, some observations from acoustic emissions indicate that shear is the dominant form of deformation, rather than extension, immediately preceding and following macroscopic failure³(p.113),¹⁹,²⁰. These analyses of acoustic emissions may detect more shear deformation than extensile deformation because the opening of extensile fractures is relatively quiet acoustically, and thus less detectable than slip-on-shear fractures. Indeed, a recent analysis highlights that acoustic emissions may not detect important kinematic deformation events during triaxial compression²¹. Moreover, comparing the numbers of acoustic emissions classified as extensile or shear may be misleading because the frequency-magnitude distribution of acoustic emissions in many experiments is similar to the Gutenberg-Richter relation for earthquakes³, and so a population

of acoustic emissions at a given time includes a wide range of event magnitudes. Consequently, analyses that count the number of acoustic emissions represent two or more acoustic emissions, which may differ by orders of magnitude in energy content, as equivalent contributions to the overall deformation.

Theoretical analyses also support the laboratory observations that indicate a prevalence of extensile deformation, at least early in loading, if not throughout loading in triaxial compression with confining stress below about 10 MPa e.g.,²². Analysis of the stresses that develop around circular holes and ellipses with the major axis aligned parallel to the σ₁ direction in an elastic medium under uniaxial loading shows that tensile stresses develop perpendicular to the σ₁ direction at the end of the cavity along the uniaxial loading axis, thus prompting fracture opening perpendicular to the σ₁ direction and growth parallel to this direction²³,²⁴. Increasing confining stress may reduce local transverse tensile stresses. Consequently, under relatively low confining stress, one may expect extensile fractures aligned parallel to the σ₁ direction to accommodate the majority of strain, and thereby dominate deformation. Under higher confining stresses, a higher proportion of shear fractures may develop. Consistent with this idea, some previous analyses found that increasing confining stress decreases the macroscopic dilatancy measured from the change in volume of the rock³, but other analyses did not observe this trend²⁵.

A more recent analysis compared the microstructures of rocks post-failure in axisymmetric extension and shortening experiments, and found that when σ₃ is tensile, the development of one or more extensile cracks oriented parallel to the σ₁ direction produces macroscopic failure¹⁷. As σ₃ approaches the compressive regime, but remains tensile, several transgranular extensile fractures coalesce into a fracture oriented 10° from the σ₁ direction. When all the principal stresses are compressive, many grain-scale extensile fractures coalesce into a shear fracture oriented from 10° to 30° from the σ₁ direction, and the orientation increases with the differential stress, σ<sub>D</sub>, at macroscopic failure¹⁷. This progression is consistent with the idea that extensile fracture development is favored under uniaxial conditions, and that increasing compressive confining stress reduces the transverse tensile stresses that develop at heterogeneities, and thereby promotes shear failure²³.

Measuring the orientation of system-spanning shear fractures post-failure in experiments indicates that their orientation approximately agrees with the predictions of Mohr-Coulomb theory. This theory determines the orientation of failure planes by maximizing the ratio of the shear stress to normal stress, i.e., the Coulomb shear stress, assuming negligible cohesion or inherent shear strength. Potential failure planes in intact rock are aligned relative to the σ₁ direction at π/4 - φ/2 under compression, where φ is the angle of internal friction²³(p. 95⁻⁹⁷). The internal friction coefficient can range from 0.5 to 1.5³(p. 2⁵), and thus the orientation of fractures under compression is expected to be about 17-32° from the σ₁ direction. Laboratory observations roughly agree with these predictions²⁶. However, recent observations of increasing fracture orientation with increasing σ₃ highlights that the internal friction coefficient may not be constant for all σ₃¹⁷. A few studies suggest that the cylindrical shape and roughly 2:1 aspect ratio of rock cores used in laboratory experiments, and the tendency for shear fractures to transverse diagonally from one end to the other due to the stress concentrations that develop at the interface between the triaxial apparatus and the rock core, may produce the apparent agreement between the measured and theoretically predicted orientation³,²⁷,²⁸.

Linear elastic models with arrays of fractures of the same length aligned parallel to the σ₁ direction suggest that the
<figure><img src="image_2.png" /><figcaption>**Fig. 1** Shear and extensile fracture development. a, b An explanation of the development of shear fractures from the linkage of extensile fractures aligned parallel to the $σ_1$ direction from Reches and Lockner¹⁸. a Fractures develop in an array-oriented about 20–30° from $σ_1$ because the lobes of higher tensile stress or lower compressive stress (light orange) develop at this orientation, and may overlap (dark orange). b However, this argument suggests that fractures could develop in quasi-parallel arrays because the lobes of higher tensile stresses produce two maxima to the left and right of the preexisting fracture tip. c, d Preexisting idealization of the evolution of deformation with increasing differential stress, $σ_D$. c At lower differential stress, most of the fractures develop quasi-parallel to $σ_1$, producing extensile-dominated deformation. d At higher differential stress, and following the peak differential stress, the macroscopic fractures are aligned near 30° from $σ_1$, indicating shear-dominated deformation.</figcaption></figure>

observed orientation of macroscopic shear fractures around 20–30° from the $σ_1$ direction arises from the tensile stress concentrations that develop at the upper and lower tips of these fractures as they open¹⁸. Because fracture opening produces a region of maximum tension, or minimum compression, inclined around 20–30° from the $σ_1$ direction at the upper and lower tips of these vertical fractures, depending on the ratio of the spacing to fracture length, the macroscopic shear fracture forms along 20–30° from the $σ_1$ direction (Fig. 1a). Reches and Lockner¹⁸ proposed that this result provides a physical interpretation of the

internal friction, the empirical parameter of the Coulomb failure criterion. However, the tensile stress concentrations that develop at the upper and lower tips of fractures aligned parallel to the $σ_1$ direction produce maximum values on both the left and right sides of the fractures. Consequently, following the criterion that fractures develop at the location of the maximum tensile stresses proposed by ref. ¹⁸, an array of parallel fractures may not develop along the expected 20–30°, but instead form subvertical arrays composed of fractures that step to the left and then to the right of the previous fracture (Fig. 1b). If one or a few of such fractures are longer than the others, then the local stress field at the tips could become asymmetric and then favor fracture development along one orientation or the other.

In summary, Griffith analyses of elliptical cavities and laboratory observations indicate that extensile deformation may dominate rock deformation at least at the onset of loading in triaxial compression for low to moderate confining stresses, if not throughout loading for confining stresses less than 10 MPa (e.g., Fig. 1c)³,¹⁷,²²⁻²⁴. Then, the Coulomb criterion and laboratory observations indicate that the macroscopic, system-spanning fracture network that is observed after the maximum stress and macroscopic failure is generally aligned between 20° and 30° from $σ_1$ for these loading conditions (e.g., Fig. 1d). Between the onset of loading and the maximum stress, different analyses indicate varying proportions of the number of tensile- and shear-dominated acoustic emissions. Analyses grounded in linear elastic fracture mechanics have used the development of mixed-mode (tensile and shear) wing cracks from the tips of inclined shear fractures to understand the stable propagation and coalescence of fractures under compression³. Although they may approximate the overall effect of a population of shear and extensile fractures, such idealized wing cracks are rarely observed in experimental data (e.g., Fig. 2)²⁹. Consequently, there is a fundamental gap in our understanding of how fractures that may nucleate as extensile fractures, develop and coalesce into system-scale fracture networks that are aligned along the orientations predicted by the Coulomb criterion.

**Observations from in situ X-ray tomography triaxial compression experiments.** Advanced experiments with in situ X-ray tomography have provided unparalleled insights into the factors that control fracture network development⁶,³⁰,³¹. Recent work has examined the evolving accumulation of the local strain in triaxial compression experiments on a wide range of rocks³². This analysis provided insight into the correlations between different types of deformation, such as the strong correlations between increases in the magnitudes of the dilative and shear strain. However, previous work has not been able to directly compare the contribution of shear and extensile deformation to the system.

To directly compare the evolving contribution of extensile and shear deformation, we examine the geometry of the fractures observed in X-ray tomography triaxial compression experiments. We perform a set of experiments with systematically varying confining stress and fluid pressure in relatively intact and heat-treated (damaged) Westerly granite. We then track the orientation of the fractures throughout loading until macroscopic failure in order to characterize the relative contribution of shear and extensile deformation (e.g., Fig. 2). Previous analyses of X-ray tomography experiments have compared the orientation of fractures with increasing differential stress using the number of fractures with a given orientation in histograms, for example (e.g.,⁴,⁵). However, this method does not provide the best representation of the fracture network because the length and volume of fractures can vary by orders of magnitude, particularly at high differential stresses close to macroscopic failure⁴⁻⁶,³³.
<figure><img src="image_3.png" /><figcaption>**Fig. 2** Fractures in X-ray tomography scans. Two-dimensional slices and three-dimensional volume renderings of fractures identified in X-ray tomography scans of experiments performed on intact (a) and damaged granite (b) with confining stress, $\sigma_2 = 15$ MPa, and fluid pressure, $P_F = 10$ MPa, at 99.9% of the differential stress at macroscopic failure. The largest fractures develop parallel to the direction of the maximum compression direction, $\sigma_1$, indicative of extensile-dominated deformation (orange arrows in a). Some fractures are aligned more obliquely to the direction of $\sigma_1$, indicative of shear-dominated deformation (blue arrow in b).</figcaption></figure>

Consequently, one very large fracture and one much smaller fracture are represented as equivalent elements of deformation in such histograms.

In the present analysis, we use the volume of the fractures when we compare the relative proportions of shear and extensile deformation. In particular, we compare the volume of fractures aligned along the orientations expected for extensile deformation, parallel and sub-parallel to the $\sigma_1$ direction, to the volume of fractures aligned along the orientations expected for shear following the Mohr-Coulomb criterion, about 17–32° from the $\sigma_1$ direction, depending on the internal friction coefficient, throughout six X-ray tomography triaxial compression experiments on intact and damaged Westerly granite cylinders with confining stresses, $\sigma_2$, of 5–20 MPa, fluid pressures, $P_F$, of zero to 10 MPa, and effective stresses, confining stress minus fluid pressure, of 5 MPa and 10 MPa. This suite of experiments provides a unique dataset that reveals differences in the deformation mechanisms of intact and damaged low-porosity crystalline rock, and rocks deformed with pore fluid at different confining stresses and the same effective stress.

## Results

**The transition from shear- to extensile-dominated deformation.** To characterize the micromechanisms of deformation in these experiments, we compare the volume of fractures that are oriented within different ranges throughout each experiment. Supplementary Figs. S1 and S2 show the axial strain-differential stress conditions when each tomogram was acquired and the total volume of fractures calculated from each tomogram. Following Mohr-Coulomb theory, using internal friction coefficients of 0.5–1.5³(p. 25), the orientation of shear fractures under compression is expected to be about 17–32°. Using the criterion for intact rocks is appropriate in this analysis because the fractures that develop preceding macroscopic failure have recently nucleated and likely do not host significant slip³⁴. Previous work suggests that the maximum internal friction coefficient may be 1.0²³, rather than 1.5, as suggested by ref. ³, producing an expected range for shear fractures of 23–32°. Consequently, we may consider fractures aligned at 17–32° or 23–32° as accommodating

shear, and fractures oriented at <17° or <23° as more vertically aligned, and thus dominated by extension rather than shear. We consider several combinations of the angle ranges used to classify each fracture as extensile- or shear-dominated. The trends that we highlight for the combination in which extensile-dominated fractures have orientations, $\theta_e = 0–17°$, and shear-dominated fractures have orientations, $\theta_s = 17–32°$, also occur for the other combinations (e.g., Supplementary Figs. S4, S5).

In the final 99% of the $\sigma_D$ at failure, all of the experiments have a larger volume of extensile fractures, $v_e^\theta$, than shear fractures, $v_s^\theta$, and thus macroscopic failure is dominated by extensile deformation rather than shear (Fig. 3). For five of the six experiments, the ratio, $r = v_e^\theta/v_s^\theta$, is generally below one until $\sigma_D$ is greater than 95–99% of the $\sigma_D$ at failure, $\sigma_D^F$, and then $r$ increases above one in the final 1–2 MPa preceding macroscopic failure. We track this ratio because the volume of both the extensile and shear fractures increase and accelerate toward failure (Figs. S3, S4). Consequently, in all but one of the experiments, shear dominates deformation until immediately preceding failure, and then extensile deformation dominates. The experiment that deviates from this trend experienced the lowest confining stress and no fluid pressure (WG10) (Fig. 3a). In this experiment, $r$ is less than one for only two tomograms, and thus extension dominates deformation throughout loading, as well as immediately preceding macroscopic failure. In all of the experiments, including WG10, the volume of extensile fractures exceeds the volume of shear fractures for the majority of the accumulated time, in terms of stress, when the normalized $\sigma_D$ is greater than 99%. In particular, when the normalized $\sigma_D$ is greater than 99%, the proportion of the differential stress in which $r$ is greater than one ranges from 75 to 100%. We characterize this dominance of the extensile deformation preceding failure in further detail in the supplementary information (Supplementary Note 1, Figs. S5, S6).

**A physical interpretation of the internal friction coefficient.** The experimental results highlight a fundamental difference between the orientation of the system-spanning fracture network that is observed following the peak stress, and the orientations of the individual fractures that propagate and coalesce with
<figure><img src="image_4.png" /><figcaption>**Fig. 3** The evolution of extensile and shear fracture development. Evolution of the ratio of the volume of the extensile-dominated ($\theta = 0-17°$) and shear-dominated ($\theta = 17-32°$) fractures, $v_\theta$, as a function of the differential stress, $\sigma_D$, normalized to the differential stress at macroscopic failure, $\sigma_D^f$, for each experiment: **a** $\sigma_2 = 5$ MPa, $P_F = 0$, intact, WG10, **b** $\sigma_2 = 10$ MPa, $P_F = 5$ MPa, intact, WG19_18, **c** $\sigma_2 = 10$ MPa, $P_F = 5$ MPa, damaged, WG06_02, **d** $\sigma_2 = 15$ MPa, $P_F = 10$ MPa, intact, WG20_19, **e** $\sigma_2 = 15$ MPa, $P_F = 10$ MPa, damaged, WG18_14, **f** $\sigma_2 = 20$ MPa, $P_F = 10$ MPa, intact, WG14. The red curve is a median filter of the data with a width of five data points. The plots on the left show the data when $\sigma_D/\sigma_D^f > 0.75$, and the plots on the right with the gray background show the data when $\sigma_D/\sigma_D^f > 0.94$. When the $v_\theta$ of the extensile fractures, $v_\theta^e$, is greater than the $v_\theta$ for the shear fractures, $v_\theta^s$, then $r = v_\theta^e/v_\theta^s$ is greater than one, and the rock hosts a larger volume of extensile fractures than shear fractures. Each pair of plots shows the evolution for a different experiment and the experimental conditions are listed at the top of each pair. The circles in the plots on the right are colored by the $\sigma_D/\sigma_D^f$, with lighter colors closer to failure.</figcaption></figure>

increasing differential stress during loading, and preceding peak stress. Previous analyses indicate that the system-spanning fracture network observed following peak stress generally follows the predictions of the Coulomb criterion. In the present experiments, the fractures that develop preceding peak stress evolve from shear- to extensile-dominated. Immediately preceding macroscopic failure, at the stress conditions that determine the Coulomb criterion, a larger volume of the fractures have orientations that are more consistent with extensile-dominated deformation than shear-dominated deformation.

Reches and Lockner¹⁸ used the geometry of lobes of higher tensile stresses that develop near the tips of fractures aligned parallel to the $\sigma_1$ direction (vertical) to provide a physical explanation of the internal friction coefficient. However, as described above, this analysis implies that an array of fractures may develop at the orientation expected from the Coulomb criterion, as well as along a subvertical arrangement, with left- and right-stepping fractures (Fig. 1a, b). More recent work has also conceptualized the development of macroscopic shear fractures from a convenient arrangement of extensile fractures

e.g.,¹⁷, but these analyses do not attempt to explain why the extensile fractures happen to align themselves along 20–30° from the $\sigma_1$ direction.

Following this work, and our experimental results, we employ linear elastic, Boundary Element Method, two-dimensional, plane strain numerical models to understand why simple fracture networks can evolve from shear- to extensile-dominated. Because we predict the geometry of fracture growth using energy optimization, these models may provide a new physical interpretation of the internal friction coefficient.

To investigate the general conditions that produce a transition from shear to tensile deformation, we model a wide range of simple initial fracture network geometries that capture the key aspects of the experimental network, and do not rely on the specific fracture arrangement observed in the experiments to produce the observations. The networks explore three key geometric parameters: the crack orientation relative to the $\sigma_1$ direction, c, the macroscopic alignment of the fractures, a, and combination of the sign of the dips of each fracture, or geometry, g, of six fractures (Fig. S7). For example, when $g = 6$, all of the
<figure><img src="image_5.png" /><figcaption>**Fig. 4** Constraints from linear elastic numerical models on the geometry of the fracture network in rocks under triaxial compression at an early stage of loading. **a** Three-dimensional plot shows the ΔW<sub>ext</sub> produced by models with selected crack orientations, **c**, macroscopic alignments, **a**, and geometries, **g**, combinations of the signs of the dips of the fractures. The fracture geometries that produce the largest ΔW<sub>ext</sub> are the most mechanically efficient, and thus most likely to form. The colors of the symbols indicate the geometries, **g**. The most efficient geometry is highlighted with red circle in the three-dimensional plot, and shown in (b). **c** It shows an example of an inefficient geometry. The blue text shows the ΔW<sub>ext</sub> for the geometries in **b**, **c**).</figcaption></figure>

fractures dip to the right, and when $g = 7$, all of the fractures dip to the left. We then calculate the change in external work, $\Delta W_{ext}$, produced by the addition and subsequent slip and opening on these fractures (Fig. 4). Under stress loading conditions, the fracture geometry that produces the largest $\Delta W_{ext}$ is the most mechanically efficient geometry, and thus the most likely to develop following the concept that fracture networks develop in order to optimize their mechanical efficiency<sup>11</sup>. Using a similar parametric approach, a previous analysis found that the geometry of thrust faults in numerical accretionary wedges that optimize work matches the geometry of these faults in laboratory experiments<sup>35</sup>. This and other previous work indicate that fault networks evolve in order to optimize the total mechanical efficiency of the system<sup>36-43</sup>.

We design the simulations to represent an idealized early stage of fracture growth when the first few fractures develop. We build rectangular numerical models with the applied effective stress of the majority of the experiments, $\sigma_2 = 5$ MPa, and then apply a $\sigma_1$ that allows the fractures to slip (Fig. S7). Increasing $\sigma_2$ to 10 MPa does not change the key results of the analysis.

The fracture geometry that produces the largest $\Delta W_{ext}$, and thus is the most efficient, shares key aspects with the experimental fracture networks (Fig. 4a). At this early stage of fracture development, the most efficient fracture network is not composed of fractures with sub-vertical orientations, but instead includes fractures oriented at $c = 35^\circ$ from the $\sigma_1$ direction. Fractures oriented at $c = 30^\circ$ in the numerical models also produce highly efficient models. This numerical result is consistent with the experiments; for all of the experiments except for experiment WG10, the majority of the six most volumetric fractures identified when $\sigma_D = 45$ MPa have orientations within the range expected for shear (Table S1). The most efficient fracture networks are not uniformly distributed throughout the system in a grid pattern, but instead are aligned along arrays of $a = 30-40^\circ$ from the $\sigma_1$ direction (Fig. 4). The agreement between the orientation of these fractures relative to the $\sigma_1$ direction and the orientation predicted by the Coulomb criterion for intact rock

with internal friction coefficients of 0.5–1.5 (17–32°)<sup>3</sup> suggests that a new physical interpretation of the internal friction coefficient may include work optimization. In particular, fracture development that optimizes the global mechanical energy is consistent with the orientation of failure planes within nominally intact rock determined from the internal friction coefficient.

To simulate the next stage of fracture network development, we model fracture propagation and interaction using the most efficient fracture geometry identified in Fig. 4 with the software GROW<sup>43</sup>. This software determines the direction of growth by optimizing work, that is, by identifying the fracture geometry that maximizes (or minimizes) $\Delta W_{ext}$ divided by the fracture area produced by fracture propagation, $\Delta A$, for stress (or displacement) loading conditions. Because the properties of the intact rock near the edges of growing fracture tips influence the direction of growth, we show the results of four combinations of shear strength, $\tau_0$, tensile strength, $t_0$, and internal friction coefficient, $\mu_i$, of the intact rock (Table S2). Three of the scenarios use the same $\tau_0$ and $t_0$ with varying internal friction, $\mu_i$, while the fourth uses higher $\tau_0$ and $t_0$ (Table S3). The rock near the tips of growing fractures is heavily damaged, and thus the strength of this rock is much lower than the bulk strength measured on intact granite samples. The Methods section and the supplementary information (Table S3) further justify the choice of these three strength parameters.

For the models with lower internal friction ($\mu_i = 0.65$ and $0.85$), the fractures coalesce and form a system-spanning fault with two branches that extend to the bottom boundaries of the models under the applied $\sigma_1$ (Fig. 5a, b). Consequently, in order to simulate additional fracture development, we identify the position and geometry of a set of four new fractures using the same approach as in Fig. 4 (Fig. S8). We refer to these different loading steps as the first and second stages of growth. For the models with $\mu_i = 0.65$ and $0.85$, the most efficient geometry identified in the second stage of growth does not include an aligned array of fractures, as in the first stage (Fig. 4), but instead includes a diffuse distribution of fractures in a grid pattern (Figs. S8, 5e, f). This result indicates that fracture networks can
<figure><img src="image_6.png" /></figure>

**Fig. 5** Numerical model results. Fracture geometries propagated using the optimization of work with different combinations of shear strength, τ₀, tensile strength, t₀, and internal friction coefficient, µᵢ, from the fracture network identified as most efficient in Fig. 4, with the parameters c=35°, a=35°, and g=6: **a-d** the first stage of growth, **e** the second stage of growth, and **i** the evolution of the ratio of the extensile/shear fracture length for each model. The models have the parameters: **a**, **e**, τ₀ = 10 MPa, t₀ = 2 MPa, µᵢ = 0.65; **b**, **f**, τ₀ = 10 MPa, t₀ = 2 MPa, µᵢ = 0.85; **c**, **g**, τ₀ = 10 MPa, t₀ = 2 MPa, µᵢ = 1.00; and **d**, **h**, τ₀ = 25 MPa, t₀ = 5 MPa, µᵢ = 1.00. Increasing the strength parameters requires increasing the σ₁ to allow growth from the tips of fractures, and so the applied σ₁ varies among the models. The red line segments show the new fractures propagated using the optimization of efficiency, and the black line segments show the initial, preexisting fracture geometry. The shaded gray areas and gray text show the orientation of fractures expected from the µᵢ using the Mohr-Coulomb criterion. **i** For all the combinations of parameters, the fracture networks develop an increasing proportion of extensile/shear fractures, similar to the experiments. The gap between the data of the first and second stage of growth of the models highlighted with red arrows reflects the four new fractures added to these models preceding the second stage of growth, following the analysis shown in Fig. S8.

favor delocalization rather than the more localized development of an aligned array of fractures conjugate to the preexisting through-going fault.

For the models with μᵢ = 1.0, the fractures propagate for several steps, but then stop growing under the applied σ₁ (Fig. 5c, d). To simulate additional fracture development, we then increase the applied σ₁ and run the models with the fracture geometries produced in the first step (Fig. 5g, h).

For all the combinations of strength parameters, the fracture networks develop an increasing proportion of extensile/shear fractures, similar to the experiments (Fig. 5i). The high temporal resolution provided by the models shows that the fracture networks experience some temporary phases in which the ratio of the extensile/shear fracture length decreases due to fracture propagation and coalescence. However, the overall trends in both the models and experiments show an increasing length or volume proportion of extensile fractures relative to shear fractures.

The models show that increasing the internal friction coefficient of the damaged, but nominally intact rock, near

growing fracture tips increases the proportion of extensile fractures that develop at the end of the simulations. Similarly, increasing the magnitude of the shear and tensile strength increases the proportion of extensile fractures. Greater shear and tensile strength of the intact rock thus favors extensile failure rather than shear failure. Moreover, the geometry of the fractures that develop in the model with the highest tested intact rock strength, τ₀ = 25 MPa, t₀ = 5 MPa, µᵢ = 1.00 (Fig. 5d, h), appears similar to the geometry of wing cracks that grow from inclined, precut fractures in experiments⁴⁴, and in numerical models⁴⁵. In contrast, lower intact rock strength promotes shear failure that links neighboring fractures in more planar geometries, and does not produce wing cracks (Fig. 5a, e).

## Discussion

Our experimental and numerical results indicate that shear is the dominant form of deformation preceding about 90–95% of the differential stress at macroscopic failure (Figs. 3, S4). After this
<figure><img src="image_7.png" /><figcaption>**Fig. 6 Summary of key results.** a Transition from shear- to extensile-dominated deformation observed in X-ray tomography triaxial compression experiments (a) and in numerical models that propagate fractures using energy optimization (b–e). a Increasing differential stress produces a larger volume of fractures oriented within the ranges expected for extensile deformation than shear deformation. b, c The most efficient numerical fracture geometries in the first stage of fracture growth do not include quasi-vertical fractures, but instead fractures oriented at 30–35° from the maximum compression direction, similar to the experiments, and consistent with the Mohr-Coulomb failure criterion for the applied internal friction coefficient. d The most efficient numerical fracture geometries in the second stage of fracture growth do not include an aligned array of fractures, but instead a distributed grid of fractures, indicating that fracture networks may delocalize. e The fracture network that develops from these distributed fractures form networks with an increasing proportion of extensile fractures, consistent with the experimental observations.</figcaption></figure>

point, extensile deformation dominates the system, in terms of the volume of fractures aligned within the ranges of orientations for extensile and shear fractures expected from Griffith analyses and the Coulomb criterion (Figs. 3, 6a). This result differs from the evolution inferred from theory and some laboratory observations of triaxial compression with low and moderate confining stresses that suggest a transition from extensile-dominated deformation under lower differential stress early in loading to shear-dominated deformation under higher differential stress later in loading, and particularly following the peak stress (Fig. 1c, d). The main difference between this previous work and the present analysis is that we are able to detect the small shear fractures early in loading, and compare the volume of the extensile and shear fractures, rather than the number of acoustic emissions, for example e.g., 19,20. Moreover, we do not observe fracture development following the peak stress, when others have inferred or directly observed the coalescence of arrays of extensile fractures into system-spanning shear fractures e.g., 18,21.

The observed evolution from shear to extensile deformation approaching failure agrees with machine learning analyses that use strain data from triaxial compression experiments that show that the local dilative strain provides more accurate predictions about the timing of macroscopic failure than the local shear strain<sup>46</sup>, and with observations of accelerating geophysical signals preceding some large earthquakes in the crust and stick-slip in the laboratory that are associated with dilatancy, such as changes in seismic wave speeds and anisotropies<sup>47,48</sup>. Microstructural observations acquired by serial thin-sectioning suggest that the first stress-induced fractures develop at preexisting intergranular and healed transgranular fractures in granite<sup>22</sup>. The approximately isotropic shapes of the most volumetrically-abundant minerals in granite (quartz and feldspar) may produce relatively random orientations of these early fractures, and not a preferred orientation parallel to the σ<sub>1</sub> direction, consistent with our observations of a higher volume proportion of shear fractures than extensile fractures early in loading (Fig. 3). Moreover, examination of the fracture network geometries at 97% of the differential stress at failure, when shear deformation dominates, show that many fractures nucleate at the edges of mineral grains, and then propagate along orientations expected of shear fractures (Fig. 7).

The evolution observed in the present analysis (Figs. 3 and 7) agrees with microstructural observations that suggest that after about 75% of the maximum differential stress, new transgranular

fractures develop at high angles to interfaces of different minerals and generally subparallel to the σ₁ direction<sup>22</sup>. Moreover, the dominance of extensile deformation immediately preceding macroscopic failure agrees with more recent observations of fractures in X-ray tomography experiments that indicate a higher number of vertically-aligned fractures toward macroscopic failure<sup>4</sup>. Here, we show that a larger volume of extensile fractures than shear fractures develops preceding macroscopic failure in both laboratory experiments and numerical models.

The increasing volume of extensile deformation toward failure agrees with the idea that the development of a through-going fracture includes the buckling and rotation of narrow columns separated by arrays of fractures aligned subparallel to the σ₁ direction<sup>12,49</sup>. Examination of the fracture networks immediately preceding failure in our X-ray tomography data reveal arrays of extensile fractures that are linked by a few shear fractures, particularly for the experiment with the lowest confining stress (Fig. 8). Applying confining stresses larger than 20 MPa may suppress the development of these columns and associated extensile deformation. For example, a previous analysis of the growth of fractures from inclined preexisting cracks using energy optimization in uniaxial and biaxial simulations show that higher confining stress promotes the propagation of shear fractures parallel to the preexisting fracture, whereas wing cracks develop under lower confining stress, and accommodate both opening and shear<sup>45</sup>.

Even if the results of the present study only apply for confining stresses lower than 20 MPa, and effective stresses lower than 10 MPa, they highlight the inadequacy of the Mohr-Coulomb criterion under these stress conditions, in agreement with recent experimental work<sup>17</sup>. This criterion derives the orientation of a potential failure plane by maximizing the Coulomb shear stress using the internal friction coefficient, an empirical value derived from the slope of the Mohr envelope. Previous analyses have argued that the internal friction coefficient should not be considered similar to the friction between two surfaces because significant sliding may only occur on fractures following macroscopic failure<sup>34</sup>. More recent work suggests that the internal friction coefficient in rocks is not only a function of friction on sliding surfaces, but rather a product of the global energy dissipation during faulting<sup>50</sup>. However, experimental and numerical analyses suggest that the frictional sliding of fractures preceding macroscopic failure contributes to the internal friction coefficient<sup>51</sup>, and thus the internal friction is some function of the
<figure><img src="image_8.png" /><figcaption>**Fig. 7** Shear and extensile fractures in X-ray tomography data. Fractures observed in two-dimensional slices of X-ray tomography data at 97% of the differential stress at failure, when the volume of shear fractures is larger than the volume of extensile fractures for experiment WG20_19 (a), and WG14 (b). The gray colors indicate the relative mass density, and thus reveal the locations of fractures and minerals. Orange and yellow arrows highlight the extensile and shear fractures, respectively. Many fractures nucleate at the edges of grains and then propagate along orientations expected for shear.</figcaption></figure>

strength of intact regions between fractures and frictional sliding on fractures⁵². Consistent with this idea, our numerical results highlight that increasing the internal friction coefficient from 0.65 to 1.0 changes the geometry of the system-spanning fracture network that develops, and in particular, produces a higher proportion of extensile/shear fracture length (Fig. 5). More generally, the present work reveals that immediately preceding macroscopic failure, near the maximum stresses that produce the Mohr envelope, the fracture network is not dominated by fractures oriented in the direction that maximizes the macroscopic Coulomb shear stress. Instead, accelerating extensile fracture development characterizes the deformation near peak stress, within 95–100% of the $\sigma_D$ at failure.

Numerical modeling provides insight into the experimental observations, and a physical interpretation of the evolution from shear to extensile-dominated deformation (Fig. 5). Constraining the geometry of first few mm-scale fractures that develop during triaxial compression using energy optimization indicates that these fractures will not be aligned sub-parallel to the $\sigma_1$ direction, but instead at 30–35° from the $\sigma_1$ direction (Fig. 4), consistent with the higher proportion of shear fractures observed early in loading in the experiments (Fig. 3). Similarly, previous analyses of planar thrust faults in numerical accretionary wedges indicate that the faults that optimize the total mechanical energy also host the highest average Coulomb stress preceding fault slip³⁵. Moreover, identifying the fracture network geometry using the optimization of the total mechanical energy produces macroscopic fracture networks at the orientation predicted by the Coulomb criterion (Fig. 5a–d, Fig. 6c). This result suggests that the internal friction coefficient is linked to the efficiency of the system, and in particular, how the orientation of a macroscopic, system-spanning fracture controls the global efficiency. This work thus provides a physical interpretation of the internal friction coefficient in rocks.

The most efficient arrangement of the new fractures that develop after a through-going fracture forms is not an aligned array, but a more diffuse distribution (Figs. 5e, f, S8). This result suggests that the fracture networks can favor delocalization, rather than the more localized development of an aligned array of fractures conjugate to the preexisting through-going fault. Analyses of the evolving spatial distribution of the high local strain

and fractures revealed in X-ray tomography triaxial compression experiments show similar temporary episodes of delocalization⁵³,⁵⁴. Such delocalization is often attributed to the presence of heterogeneities⁵⁵. Several large earthquakes in southern and Baja California were preceded by the localization of low-magnitude seismicity that included phases of delocalization⁵⁶. Consequently, understanding the factors that control episodes of delocalization may provide key insights into the precursory processes that accelerate preceding large earthquakes. The present analysis indicates that delocalization can occur as fracture networks develop in order to optimize the total mechanical efficiency of the system.

The numerical models suggest that first the system develops a larger proportion of shear fractures, and then it develops a larger and increasing proportion of extensile fractures (Fig. 5i), consistent with the experimental observations (Fig. 3). Because these models simulate fracture growth and coalescence with work optimization, the results suggest that the experimentally observed transition from shear to extension occurs because the fracture network develops in order to optimize the mechanical efficiency of the system.

## Methods

**Experimental design.** We performed six triaxial compression experiments at beamline ID19 at the European Synchrotron Radiation Facility, in Grenoble, France. For each experiment, we inserted one 10 mm tall and 4 mm diameter cylinder of Westerly granite in the Hades triaxial compression apparatus⁵⁷ installed on the beamline. We imposed a confining stress (5–20 MPa) using pressurized oil against the jacket surrounding the core, included a fluid pressure of zero to 10 MPa, and increased the axial stress, $\sigma_1$, in steps of 0.5–5 MPa, with smaller steps closer to failure, until the rock failed with a stress drop (Table S4, Fig. S1). We varied the confining stress and fluid pressure so that five of the six experiments experienced the same effective stress (confining pressure minus pore fluid pressure was equal to 5 MPa). After each increase in $\sigma_1$, we acquired an X-ray scan within 1.5 min while the rock was under stress.

We deformed both intact and heat-treated (damaged) Westerly granite in order to assess the influence of preexisting damage on
<figure><img src="image_9.png" /><figcaption>**Fig. 8** Shear and extensile fractures in X-ray tomography data. Fracture networks observed in the final tomogram acquired preceding macroscopic failure, when the differential stress is within 0.1 MPa of the failure stress for experiment (a) WG20_19, (b) WG14, and (c) WG10. The orange and yellow arrows highlight the extensile and shear fractures, respectively. Fractures often initiate at the edges of minerals. Arrays of extensile fractures develop in the experiment with the lowest confining stress (c). Shear fractures extend from the tips of extensile fractures, and link neighboring extensile fractures.</figcaption></figure>

fracture development. We damaged the cores by heating them in an oven (Table S4). Westerly granite is a low-porosity crystalline rock dominated by interlocking quartz, feldspar, and biotite. The initial porosity is lower than 1% for the intact rocks. For the experiments that include fluid pressure, we saturated the granite cores in

deionized water in a vacuum chamber for two weeks preceding the experiment.

Following each experiment, we reconstruct the radiographs into three-dimensional 1600 × 1600 × 1600 voxel volumes, in which each voxel side length is 6.5 µm. During reconstruction, we apply algorithms to remove noise, such as ring artefacts. We remove remaining noise in the three-dimensional data using the software Avizo3D™, including the application of a non-local means filter⁵⁸. We then segment the solid rock from the pores and fractures using an algorithm similar to Otsu's thresholding technique to identify a global threshold between the solid material and the fractures and pores⁵⁹. This thresholding technique is not strongly influenced by noise⁵⁹ and produces segmented scans of sandstone with porosity similar to values measured with sample imbibition⁶⁰. We test the influence of varying this threshold in each experiment, and find that varying the threshold does not change the central results (Fig. S9).

We then apply several processing techniques to the segmented scans in Avizo3D™, including Label Analysis, that identify individual fractures using connected clusters of voxels, and calculate the geometric characteristics and orientation of fractures relative to the vertical and $σ_1$ direction. In particular, we use the eigenvectors of the best-fit ellipsoids of the fractures to identify the orientations. To further remove noise from the segmented fracture networks, we only analyze fractures with volumes greater than 20 voxels. In general, the fractures are much larger than this limit: the mean fracture volume in each experiment is 196, 210, and 320 voxels when the normalized differential stress is 0–50%, 50–75%, and >75%. The volume of the fractures is characterized by the number of voxels that comprise each identified, independent connected cluster of voxels. Because we use the volume of the fractures to compare the dominance of shear and extensile deformation, the most volumetric fractures, with volumes generally >200 voxels, dominantly control the resulting volume proportion of shear and extensile fractures. To remove artefacts that arise from phase contrasts at mineral boundaries, we remove identified pores with shape anisotropies (one minus the shortest/longest eigenvalue) below 0.95.

Because fractures must have apertures greater than the spatial resolution of the tomograms (6.5 µm/voxel) in order to be identified, one may suspect that an observational bias could exist that favors the detection of extensile-dominated fractures rather than shear-dominated fractures. For a perfectly planar fracture that hosts shear, no dilation is required for shear displacement. However, because real fractures have roughness, shear deformation requires some dilation in order to move the rough surfaces past each other. This shear movement may also involve the breakage of asperities that reduce the amount of dilation expected from the roughness. To estimate the extent of this bias in the tomograms, we compare the shape anisotropy (one minus the aperture divided by the length) of fractures classified as shear and extensile using one of the combinations of orientations described above. The apertures of the shear and extensile fractures are similar at a given stress step in each experiment, so we may compare the anisotropies of the fractures to assess the strength of the bias. If the observational bias was strong, we would expect larger anisotropies for the shear fractures than the extensile fractures. In contrast, we find that the average shape anisotropy of the extensile-dominated fractures ($\theta$=0–17°) is similar or slightly larger than the average anisotropy of the shear-dominated fractures ($\theta$=17–32°) throughout loading in each experiment (Fig. S10). This result indicates that this observational bias does not influence the comparison of the volume proportion of extensile and shear fractures. In addition, the high average anisotropy of the fractures suggests that the orientation of the fractures is representative of the overall orientation, and not
biased by highly tortuous fractures. Moreover, the similar shape anisotropy of these two types of fractures suggests that the volume of the fractures is the most appropriate proxy for the influence of each type of fracture on the deformation field.

As the differential stress increases, the fractures grow in length, and new fractures nucleate between preexisting fractures, thereby decreasing the average distance between fractures. If the fracture density becomes large enough, the directions of the principal stresses near a fracture may rotate away from the imposed global directions of the principal stresses. This rotation could cause some fractures that are aligned at 30° from the global σ₁ direction to be aligned at 5° from the local σ₁ direction, for example. Although these fractures may accommodate a significant proportion of opening, such fractures would then be labeled as shear fractures in the present analysis. Previous analyses suggest that when the ratio of the spacing between fractures and fracture length becomes less than one, the growth of one fracture may influence another via the perturbation of the local stress field caused by fracture growth and opening⁶¹. This perturbation may promote a rotation of the local principal stresses away from the global principal stress directions. However, X-ray tomography triaxial compression experiments show that the lengths of the fractures at one stress step vary by orders of magnitude³³, and accordingly the distances between fractures vary. Consequently, we cannot identify the time in each experiment when the fracture density is large enough to produce a significant rotation of the local principal stress directions. Moreover, the time when this occurs may differ for each fracture, depending on its length and distance to other fractures. Consequently, we follow the approach of decades of previous work that have used the orientation of fractures relative to the global directions of the principal stresses to compare the relative proportions of extensile and shear deformation e.g.,²².

We note that some fractures may be best classified as mixed-mode fractures, rather than pure mode I or mode II/III. Indeed, in order for a rough fracture to shear, some dilation may be required. Consequently, we emphasize that the fractures classified as extensile or shear using their orientation are only dominated by extensile or shear deformation, and do not exclusively host one of these endmembers of deformation.

We note that this technique cannot detect fractures when they have apertures below the spatial resolution. However, because the analysis uses the volume of fractures, rather than counting the number of fractures, the most volumetric fractures dominantly influence the results.

**Numerical model design.** To identify the arrangement of fractures in the first stage of fracture growth, we use the software Fric2D⁶². Fric2D uses continuum mechanics to calculate the stresses and strains within a homogeneous linear elastic material containing fractures and frictional faults in two-dimensional plane strain⁶². The two-dimensional plane strain models are one meter in the direction out of the board/screen, or z-direction. The x-direction is parallel to the σ₂ direction and y-direction is parallel to the σ₁ direction.

We select the loading conditions of the first set of numerical models to match the effective stress of the majority of the experiments (σ₂ = 5 MPa), and to produce slip along the fractures with the lowest amount of loading (e.g., σ₁ = 50 MPa for the lowest strength parameters of the intact rock). Models with higher applied confining stress (σ₂ = 10 MPa) produce similar results to models with σ₂ = 5 MPa. Applying the lowest σ₁ to produce failure ensures that the numerical models capture the response of the fractures to loading in a realistic manner. In an experiment, fractures begin to grow as soon as the σ₁ is sufficient

to propagate new fractures, and so applying a σ₁ larger than sufficient to produce failure is an unrealistic loading condition that will likely produce unrealistic fault geometries.

We thus prescribe loading conditions that represent an early stage of the experiments. We select the length of the numerical fractures (0.3 mm) using the lengths of the largest six fractures identified in each experiment when σD = 45 MPa (Fig. S7). We select an element length that is large enough to satisfy the requirement of the linear elastic models that slip or opening on an element must be less than half of the element length. We use a rectangular model geometry with a side length of 8 mm perpendicular to the σ₁ direction (horizontal) and 4 mm parallel to the σ₁ direction (vertical) to ensure that the system shape does not promote the development of fractures at 20–30° from the σ₁ direction, as suggested by previous work²⁷,²⁸.

To simulate the growth of fractures, we use the software GROW⁴³, which calls Fric2D. GROW accurately simulates the propagation path and linkage of two nearby fracture tips separated by a releasing step⁶³, the propagation of wing cracks from inclined fractures observed in uniaxial and biaxial compression experiments⁶⁴, the coalescence of hundreds of fractures observed in uniaxial compression experiments⁶⁴, and the development of thrust faults in accretionary prisms⁶⁵.

Recent work describes in detail the key algorithms of GROW⁶⁴,⁶⁵. In summary, in the single run mode of fracture growth, the fracture tip that grows during a given time step is the tip with the highest sum of the absolute value of the mode-I and mode-II stress intensity factors. The models simulate fracture growth through intact rock by adding elements to the tips of the preexisting faults that have the mechanical properties of the intact and damaged rock near fracture tips, including the tensile strength and shear strength. Tables S2 and S3 list the mechanical and frictional properties applied in the models, which are representative of intact, but damaged, granite and faults in granite. After GROW identifies the most efficient orientation of growth from a preexisting fracture, this intact rock element gains the mechanical and frictional properties (e.g., sliding friction) of the preexisting fractures, which matches the properties of faults in granite in these simulations. If none of the intact rock elements added to the tips of the preexisting faults fail in tension (when the tensile stress is greater than the applied tensile strength) or in shear following the Coulomb criterion, then this fracture tip does not propagate in the simulation, and stops growing. A GROW simulation ends when all of the fractures stop growing, or all of the tips intersect the boundaries of the system or other fractures.

We test different combinations of the shear strength, tensile strength, and internal friction coefficient of the intact rock near the edges of growing fractures (Table S3). Because the values determine whether the new elements added to the tips of growing fracture fail, in the models with higher strength parameters we must apply a higher σ₁ to produce fracture growth than the models with lower strength parameters, with σ₁ = 50 MPa. The values of these properties for the bulk rock that is greater than one element length away from the growing fracture tip do not influence the predicted fracture geometries.

The properties of these intact rock elements represent the nominally intact, but heavily damaged, rock at the tips of preexisting fractures. The tested internal friction coefficients are within the lower range measured for nominally intact granite⁶⁶. The minimum tested internal friction coefficient is the sliding friction coefficient of the preexisting fractures in the models (0.65) because rocks should not have an internal friction coefficient that is lower than the sliding friction coefficient. The ratio of the tensile and shear strength for low porosity, crystalline rocks is generally 1:5³. For low porosity, nominally intact, crystalline rocks such as granite, the inherent shear strength
ranges from 50 to 100 MPa, and the tensile strength ranges from about 5–20 MPa³,⁶⁷,⁶⁸. Following the approach of Fattoruso et al.⁶⁴, we selected the values of 2 MPa and 10 MPa, and 5 MPa and 25 MPa for the tensile and shear strength because the rock near the tips of growing fractures is heavily damaged, and thus the strength of this rock is much lower than the bulk tensile and shear strength measured for intact granite samples. These values determine whether a newly added element will fail at a given level of differential stress, and thus may influence the direction of growth from a propagating fracture tip.

Testing the influence of $\tau_0$, for constant $\mu_i$ and $t_0$, the run mode of GROW, and the inclusion of ten fractures at the beginning of the simulations shows that the proportion of extensile to shear fracture length increases for all of the tested parameters (Fig. S11), consistent with the experiments. The frictional and strength properties of the fractures in the models are the same for all of the fractures.

To compare the experiments and the numerical models, we calculate the total fracture length by summing the length of all the linear elements of all the fractures at a given model time step. We calculate the fracture orientation using individual elements, and so one connected fracture may have different orientations for each element. Ideally, we would have calculated the fracture orientations in the experiments with this method also, but this is not feasible. However, the shape anisotropy of the fractures in the experiments is very high (with values generally >0.95), and so there are not many highly tortuous fractures that would produce varying orientations if we measured the orientation segment by segment.

**Model limitations.** The numerical models are necessarily a simplification of the fracture development that occurs in the experiments. The two-dimensional system prevents the development of polymodal faulting⁶⁹. The low strain simplification of the linear elastic models prevents simulating the full loading of a triaxial compression experiment, and thus requires simulating snapshots of the experiment at particular stress steps. The shape of the rock cores used in the experiments produce a heterogeneous stress field that can produce stress concentrations at the upper and lower corners of the rock cores, where the rock is in contact with the triaxial deformation apparatus. Depending on the lubrication between the rock and the apparatus, these stress concentrations can promote the development of fractures from one corner to another, in an array aligned at 20–30° from the $\sigma_1$ direction²⁷,²⁸. In order to eliminate such boundary effects, we carefully designed the shape of the numerical models and the placement of the fractures in the models. Consequently, the numerical models may be considered to be a central portion of the rock core that is not strongly influenced by boundary effects. Despite these limitations, the models produce key characteristics observed in the experiments, including the preference for shear fractures early in loading, and the increasing proportion of extensile fractures relative to shear fractures. The numerical models do not include the heterogeneity that arises from the interlocking crystalline mineral structure of the granite. However, they capture key aspects of the fracture development observed in the experiments. Consequently, the microstructure of this granite may not exert a dominant control on the fundamental characteristics of fracture development observed in the models, such as the transition from shear deformation to extensile. The microstructure influences where the first fractures nucleate, but this nucleation stage is not captured in the numerical models. Rather, the earliest stage of the experiments that the models capture is when the differential stress is 45 MPa and when the longest fractures are about 0.3 mm.

## Data availability

The denoised, filtered, and segmented tomograms are available on Norstore with DOI 10.11582/2023.00007 (https://archive.norstore.no/pages/public/datasetDetail.jsf?id=10.11582/2023.00036)⁷¹. The software GROW and Fric2D are available on GitHub (https://github.com/jmcbbeck/GROW).

## Code availability

The scripts used to process the tomograms are available on Norstore with DOI 10.11582/2023.00036 (https://archive.sigma2.no/pages/public/datasetDetail.jsf?id=10.11582/2023.00036)⁷¹. The software GROW and Fric2D are available on GitHub (https://github.com/jmcbbeck/GROW).

Received: 6 March 2023; Accepted: 25 September 2023;

Published online: 05 October 2023

## References

1. Kemeny, J. M. & Cook, N. G. Micromechanics of deformation in rocks. *Toughening Mechanisms in Quasi-Brittle Materials* (pp. 155–188. Springer, Dordrecht, 1991).

2. Brace, W. F. & Bombolakis, E. G. A note on brittle crack growth in compression. *J. Geophys. Res.* **68**, 3709–3713 (1963).

3. Paterson, M. S. & Wong, T. F. *Experimental Rock Deformation: The Brittle Field* 348 (Springer, Berlin, 2005).

4. Cartwright-Taylor, A. et al. Catastrophic failure: how and when? Insights from 4-D in situ X-ray microtomography. *J. Geophys. Res. Solid Earth* **125**, e2020JB019642 (2020).

5. Renard, F. et al. Critical evolution of damage toward system-size failure in crystalline rock. *J. Geophys. Res. Solid Earth* **123**, 1969–1986 (2018).

6. Renard, F. et al. Volumetric and shear processes in crystalline rock approaching faulting. *Proc. Natl Acad. Sci.* **116**, 16234–16239 (2019).

7. McBeck, J., Kandula, N., Aiken, J. M., Cordonnier, B. & Renard, F. Isolating the factors that govern fracture development in rocks throughout dynamic in situ X-ray tomography experiments. *Geophys. Res. Lett.* **46**, 11127–11135 (2019).

8. McBeck, J. A., Aiken, J. M., Mathiesen, J., Ben-Zion, Y. & Renard, F. Deformation precursors to catastrophic failure in rocks. *Geophys. Res. Lett.* **47**, e2020GL090255 (2020).

9. Orowan, E. Fracture and strength of solids. *Rep. Prog. Phys.* **12**, 185–232 (1949).

10. Pollard, D. D. & Aydin, A. Progress in understanding jointing over the past century. *Geol. Soc. Am. Bull.* **100**, 1181–1204 (1988).

11. Cooke, M. L. & Madden, E. H. Is the Earth lazy? A review of work minimization in fault evolution. *J. Struct. Geol.* **66**, 334–346 (2014).

12. Peng, S., & Johnson, A. M. Crack growth and faulting in cylindrical specimens of Chelmsford granite. Int. J. Rock Mech. Mining Sci. Geomech. Abstr. **9**, 37–86 (1972).

13. Ashby, M. F. & Sammis, C. G. The damage mechanics of brittle solids in compression. *Pure Appl. Geophys.* **133**, 489–521 (1990).

14. Horii, H. & Nemat-Nasser, S. Compression-induced microcrack growth in brittle solids: axial splitting and shear failure. *J. Geophys. Res. Solid Earth* **90**, 3105–3125 (1985).

15. Wibberley, C. A., Petit, J. P. & Rives, T. Micromechanics of shear rupture and the control of normal stress. *J. Struct. Geol.* **22**, 411–427 (2000).

16. Heard, H. C. Transition from brittle fracture to ductile flow in Solenhofen limestone as a function of temperature, confining pressure, and interstitial fluid pressure. In: (eds. Griggs, D. T. & Handin, J.) *Rock Deformation* **79**, 193–226 (Geological Society of America Memoirs, 1960).

17. McCormick, C. A. & Rutter, E. H. An experimental study of the transition from tensile failure to shear failure in Carrara marble and Solnhofen limestone: Does "hybrid failure" exist? *Tectonophysics* **844**, 229623 (2022).

18. Reches, Z. E. & Lockner, D. A. Nucleation and growth of faults in brittle rocks. *J. Geophys. Res. Solid Earth* **99**, 18159–18173 (1994).

19. Stanchits, S., Vinciguerra, S. & Dresen, G. Ultrasonic velocities, acoustic emission characteristics and crack damage of basalt and granite. *Pure Appl. Geophys.* **163**, 975–994 (2006).

20. Graham, C. C., Stanchits, S., Main, I. G. & Dresen, G. Comparison of polarity and moment tensor inversion methods for source analysis of acoustic emission data. *Int. J. Rock Mech. Mining Sci.* **47**, 161–169 (2010).

21. Cartwright-Taylor, A. et al. Seismic events miss important kinematically governed grain-scale mechanisms during shear failure of porous rock. *Nat. Commun.* **13**, 1–14 (2022).

22. Tapponier, P. & Brace, W. F. Development of stress-induced micro-cracks in Westerly granite. *Int. J. Rock Mech. Mining Sci.* **13**, 103–112 (1976).
23. Jaeger, J. C. & Cook, N. G. W. *Fundamentals of Rock Mechanics*. 3rd edn (Chapman & Hall, London, 1979).

24. Davis, T., Healy, D., Bubeck, A. & Walker, R. Stress concentrations around voids in three dimensions: the roots of failure. *J. Struct. Geol.* **102**, 193–207 (2017).

25. Brace, W. F., Paulding, B. W. Jr & Scholz, C. H. Dilatancy in the fracture of crystalline rocks. *J. Geophys. Res.* **71**, 3939–3953 (1966).

26. Mogi, K. Some precise measurements of fracture strength of rocks under uniform compressive stress. *Felsmech. Ingenieurgeol.* **4**, 41 (1966).

27. Paul, B., & Gangal, M.. Initial and subsequent fracture curves for biaxial compression of brittle materials. In *Proc. 8th US Symposium on Rock Mechanics (USRMS)* (1966).

28. Paul, B. Macroscopic criteria for plastic flow and brittle fracture. *Fracture 2*, 313–496 (1968).

29. McBeck, J., Aiken, J. M., Cordonnier, B., Ben-Zion, Y. & Renard, F. Predicting fracture network development in crystalline rocks. *Pure Appl. Geophys.* **179**, 275–299 (2022).

30. Huang, L. et al. Synchrotron X-ray imaging in 4D: multiscale failure and compaction localization in triaxially compressed porous limestone. *Earth Planet. Sci. Lett.* **528**, 115831 (2019).

31. Liu, S. & Huang, Z. Analysis of strength property and pore characteristics of Taihang limestone using X-ray computed tomography at high temperatures. *Sci. Rep.* **11**, 1–14 (2021).

32. McBeck, J., Ben-Zion, Y. & Renard, F. The mixology of precursory strain partitioning approaching brittle failure in rocks. *Geophys. J. Int.* **221**, 1856–1872 (2020).

33. Kandula, N. et al. Dynamics of microscale precursors during brittle compressive failure in Carrara marble. *J. Geophys. Res. Solid Earth* **124**, 6121–6139 (2019).

34. Handin, J. On the Coulomb–Mohr failure criterion. *J. Geophys. Res.* **74**, 5343–5348 (1969).

35. McBeck, J. A., Cooke, M. L., Herbert, J. W., Maillot, B. & Souloumiac, P. Work optimization predicts accretionary faulting: an integration of physical and numerical experiments. *J. Geophys. Res. Solid Earth* **122**, 7485–7505 (2017).

36. Lachenbruch, A. H. & Thompson, G. A. Oceanic ridges and transform faults: their intersection angles and resistance to plate motion. *Earth Planet. Sci. Lett.* **15**, 116–122 (1972).

37. Melosh, H. J. & Williams, C. A. Jr Mechanics of graben formation in crustal rocks: a finite element analysis. *J. Geophys. Res. Solid Earth* **94**, 13961–13973 (1989).

38. Du, Y. & Aydin, A. The maximum distortional strain energy density criterion for shear fracture propagation with applications to the growth paths of en-echelon faults. *Geophys. Res. Lett.* **20**, 1091–1094 (1993).

39. Okubo, C. H. & Schultz, R. A. Evolution of damage zone geometry and intensity in porous sandstone: insight gained from strain energy density. *J. Geol. Soc.* **162**, 939–949 (2005).

40. Del Castello, M. & Cooke, M. L. Underthrusting-accretion cycle: work budget as revealed by the boundary element method. *J. Geophys. Res. Solid Earth* **112**, 1–14 (2007). B12404.

41. Olson, E. L. & Cooke, M. L. Application of three fault growth criteria to the Puente Hills thrust system, Los Angeles, California, USA. *J. Struct. Geol.* **27**, 1765–1777 (2005).

42. Griffith, W. A. & Cooke, M. L. Mechanical validation of the three-dimensional intersection geometry between the Puente Hills blind-thrust system and the Whittier fault, Los Angeles, California. *Bull. Seismol. Soc. Am.* **94**, 493–505 (2004).

43. McBeck, J. A., Madden, E. H. & Cooke, M. L. Growth by Optimization of Work (GROW): a new modeling tool that predicts fault growth through work minimization. *Comput. Geosci.* **88**, 142–151 (2016).

44. Bobet, A. & Einstein, H. H. Fracture coalescence in rock-type materials under uniaxial and biaxial compression. *Int. J. Rock Mech. Mining Sci.* **35**, 863–888 (1998).

45. Madden, E. H., Cooke, M. L. & McBeck, J. Energy budget and propagation of faults via shearing and opening using work optimization. *J. Geophys. Res. Solid Earth* **122**, 6757–6772 (2017).

46. McBeck, J., Aiken, J. M., Ben-Zion, Y. & Renard, F. Predicting the proximity to macroscopic failure using local strain populations from dynamic in situ X-ray tomography triaxial compression experiments on rocks. *Earth Planet. Sci. Lett.* **543**, 116344 (2020).

47. Shreedharan, S., Bolton, D. C., Rivière, J. & Marone, C. Competition between preslip and deviatoric stress modulates precursors for laboratory earthquakes. *Earth Planet. Sci. Lett.* **553**, 116623 (2021).

48. Chiarabba, C., De Gori, P., Segou, M. & Cattaneo, M. Seismic velocity precursors to the 2016 Mw 6.5 Norcia (Italy) earthquake. *Geology* **48**, 924–928 (2020).

49. Wong, T. F. Micromechanics of faulting in Westerly granite. *Int. J. Rock Mech. Mining Sci. Geomech. Abstr.* **19**, 49–64 (1982).

50. Reches, Z. E. & Wetzler, N. An energy-based theory of rock faulting. *Earth Planet. Sci. Lett.* **597**, 117818 (2022).

51. Schoepfer, M. P., Childs, C. & Manzocchi, T. Three-dimensional failure envelopes and the brittle-ductile transition. *J. Geophys. Res. Solid Earth* **118**, 1378–1392 (2013).

52. Savage, J. C., Byerlee, J. D. & Lockner, D. A. Is internal friction friction? *Geophys. Res. Lett.* **23**, 487–490 (1996).

53. McBeck, J., Ben-Zion, Y. & Renard, F. Volumetric and shear strain localization throughout triaxial compression experiments on rocks. *Tectonophysics* **822**, 229181 (2022).

54. McBeck, J., Ben-Zion, Y. & Renard, F. Fracture network localization preceding catastrophic failure in triaxial compression experiments on rocks. *Front. Earth Sci.* **9**, 778811 (2021).

55. Alava, M. J., Nukala, P. K. & Zapperi, S. Statistical models of fracture. *Adv. Phys.* **55**, 349–476 (2006).

56. Ben-Zion, Y. & Zaliapin, I. Localization and coalescence of seismicity before large earthquakes. *Geophys. J. Int.* **223**, 561–583 (2020).

57. Renard, F. et al. A deformation rig for synchrotron microtomography studies of geomaterials under conditions down to 10 km depth in the Earth. *J. Synchrotron Radiat.* **23**, 1030–1034 (2016).

58. Buades, A., Coll, B. & Morel, J. M. A non-local algorithm for image denoising. In *Proc. IEEE Computer Society Conference on Computer Vision and Pattern Recognition (CVPR’05)*, 2, 60–65 (IEEE, 2005)..

59. McBeck, J. A., Cordonnier, B. & Renard, F. The influence of spatial resolution and noise on fracture network properties calculated from X-ray microtomography data. *Int. J. Rock Mech. Mining Sci.* **147**, 104922 (2021).

60. Renard, F. et al. Dynamic in situ three-dimensional imaging and digital volume correlation analysis to quantify strain localization and fracture coalescence in sandstone. *Pure Appl. Geophys.* **176**, 1083–1115 (2019).

61. Olson, J. E. Joint pattern development: effects of subcritical crack growth and mechanical crack interaction. *J. Geophys. Res. Solid Earth* **98**, 12251–12265 (1993).

62. Cooke, M. L. & Pollard, D. D. Bedding-plane slip in initial stages of fault-related folding. *J. Struct. Geol.* **19**, 567–581 (1997).

63. McBeck, J., Cooke, M. & Madden, E. Work optimization predicts the evolution of extensional step over steps within anisotropic host rock: Implications for the San Pablo Bay, CA. *Tectonics* **36**, 2630–2646 (2017).

64. Fattaruso, L., Cooke, M. & McBeck, J. The influence of fracture growth and coalescence on the energy budget leading to failure. *Front. Earth Sci.* **10**, 1–19 (2022). 853030.

65. McBeck, J., Cooke, M. & Fattaruso, L. Predicting the propagation and interaction of frontal accretionary thrust faults with work optimization. *Tectonophysics* **786**, 228461 (2020).

66. Scholz, C. H.. The Mechanics of Earthquakes and Faulting. (Cambridge University Press, 2019)

67. Zhang, Z. X. An empirical relation between mode I fracture toughness and the tensile strength of rock. *Int. J. Rock Mech. Mining Sci.* **39**, 401–406 (2002).

68. Liu, Y. et al. Tensile strength and fracture surface morphology of granite under confined direct tension test. *Rock Mech. Rock Eng.* **54**, 4755–4769 (2021).

69. Healy, D. et al. Polymodal faulting: time for a new angle on shear failure. *J. Struct. Geol.* **80**, 57–71 (2015).

70. McBeck, J. A. Granite experiments [Data set]. Norstore. https://doi.org/10.11582/2023.00007 (2023).

71. McBeck, J. A. Processing scripts for tomography data [Data set]. Norstore. https://doi.org/10.11582/2023.00036 (2023).

Acknowledgements

This study was funded by the Research Council of Norway (project no. 300435 to J.M.), UNINETT Sigma2 AS (project no. NN9806K to J.M. and project no. NS9073K to F.R.), and the European Research Council (project no. 101019628 to F.R.). Paul Meakin provided useful suggestions on this manuscript. We thank Editor Joe Aslin, Alexis Cartwright-Taylor, and Christoph Schrank for thorough reviews that improved this manuscript.

Author contributions

J.M. designed the experimental and numerical analysis and wrote the initial drafts of the manuscript. L.F. and M.C. contributed numerical modeling assistance and recommendations on the scope of the analyses, and edited the manuscript. B.C. helped perform the experiments and edited the manuscript. F.R. helped perform the experiments, guide the scope of the analysis, and edited the manuscript.

Competing interests

The authors declare no competing interests.
**Additional information**

**Supplementary information** The online version contains supplementary material available at https://doi.org/10.1038/s43247-023-01023-w.

**Correspondence** and requests for materials should be addressed to Jess McBeck.

**Peer review information** *Communications Earth & Environment* thanks Alexis Cartwright-Taylor, and Christoph Schrank for their contribution to the peer review of this work. Primary Handling Editor: Joe Aslin. A peer review file is available.

**Reprints and permission information** is available at http://www.nature.com/reprints

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

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

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2023