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

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


Introduction
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 [1][2][3].Recent work suggests that fracture network development controls the predictability of failure [4][5][6].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 [7][8].Theories from linear elastic fracture mechanics predict with success the conditions for propagation, and geometry of isolated, individual fractures [9][10].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 compression [e.g., 2,11].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 heterogeneous 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 [11].1.1.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 [3].We adopt here the terminology of "extensile fractures" because "extension" refers to strain, whereas "tension" refers to stress [e.g., 3].
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 [3].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 [3].Laboratory analyses suggest that the systemspanning fracture networks that accommodate shear develop from the coalescence of preexisting extensile fractures [12][13][14][15][16][17], and this process may include the buckling of columns between the extensile fractures [12].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 [18].
However, some observations from acoustic emissions indicate that shear is the dominant form of deformation, rather than extension, at higher differential stresses and thus immediately preceding macroscopic failure [ref.3 p.113, 19-20].This apparent difference in the dominance of shear and extensile failure may occur 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 [21].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 [3], 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., 22].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 [23,24].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 [3], but other analyses did not observe this trend [25].
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 [17].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-30° from the  !direction, and the orientation increases with the differential stress at macroscopic failure [17].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 [23].
Measuring The internal friction coefficient can range from 0.5-1.5 [3 p. 25], 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 [26].However, recent observations of increasing fracture orientation with increasing  " highlights that the internal friction coefficient may not be constant for all  " [17].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 [3,[27][28].
Linear elastic models with arrays of fractures of the same length aligned parallel to the  !direction suggest that the observed orientation of macroscopic shear fractures around 20-30° from the  !direction arises from the tensile stress concentrations that develop at the upper and lower tips of these fractures as they open [18].Because fracture opening produces a region of maximum tension, or minimum compression, inclined around 20-30° from the  !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  !direction (Figure 1a).Reches & Lockner [18] 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  !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 Reches & Lockner [18], 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 (Figure 1b).
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, if not throughout loading for confining stresses less than 10 MPa (e.g.,

Figure 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-30° from  ! for these loading conditions (e.g., Figure 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 [3].However, such idealized wing cracks are rarely observed in experimental data (e.g., Figure 2) [e.g., 29].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 [6,[30][31].Recent work has examined the evolving accumulation and partitioning of the local strain tensors in triaxial compression experiments on a wide range of rocks [32].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., Figure 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., [4][5].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 [4][5][6]33].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 explicitly 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,

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.
The supplementary Figure S1 and Figure 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 [ref.3 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 [34].Previous work suggests that the maximum internal friction coefficient may be 1.0 [23], rather than 1.5, as suggested by [3], producing an expected  In the final 99% of the  ' at failure, all of the experiments have a larger volume of extensile fractures,  + ) , than shear fractures,  + * , and thus macroscopic failure is dominated by extensile deformation rather than shear (Figure 3).For five of the six experiments, the ratio,  =  + ) / + * , is generally below one until  ' is greater than 95-99% of the  ' at failure,  ' ( , and then  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 (Figure S3, Figure 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) (Figure 3a).In this experiment,  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  ' is greater than 99%.In particular, when the normalized  ' is greater than 99%, the proportion of the differential stress in which  is greater than one ranges from 75-100% of the total differential stress for each experiment.We characterize this dominance of the extensile deformation preceding failure in further detail in the supplementary information (Text S1, Figure S5, Figure 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 increasing differential stress during loading.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 & Lockner [18] used the geometry of lobes of higher tensile stresses that develop near the tips of fractures aligned parallel to the  !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 (Figure 1a-b).More recent work has also conceptualized the development of macroscopic shear fractures from a convenient arrangement of extensile fractures [e.g., 17], but these analyses do not attempt to explain why the extensile fractures happen to align themselves along 20-30° from the  !direction.
Following this work, and our experimental results, we employ linear elastic, Boundary Element Method, two-dimensional, plane strain numerical models to understand why the fracture networks 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 simulate fracture growth, we consider a wide range of initial fracture network geometries using three geometric parameters: the crack orientation relative to the  !direction, , the macroscopic alignment of the fractures, , and combination of the sign of the dips of each fracture, or geometry, , of six fractures (Figure S7).We then calculate the change in external work, ∆ ),-, produced by the addition and subsequent slip and opening on these fractures (Figure 4).Under stress loading conditions, the fracture geometry that produces the largest ∆ ),-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 [11].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 [35].This and other previous work indicate that fault networks evolve in order to optimize the total mechanical efficiency of the system [36][37][38][39][40][41][42][43].
We design the simulations to represent an 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,  & = 5 MPa, and then apply a  ! that allows the fractures to slip (Figure S7).Increasing  & to 10 MPa does not change the key results of the analysis.
The fracture geometry that produces the largest ∆ ),-, and thus is the most efficient, shares key aspects with the experimental fracture networks (Figure 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 =35° from the  !direction, consistent with the experiments.The most efficient fracture networks are not uniformly distributed throughout the system in a grid pattern, but instead are aligned along arrays of =30-40° from the  !direction (Figure 4).The agreement between the orientation of these fractures relative to the σ !direction and the orientation predicted by the Coulomb criterion for intact rock with internal friction coefficients of 0.5-1.5 (17-32°) [3] 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 Figure 4 with the software GROW [43].This software determines the direction of growth by optimizing work, that is, by identifying the fracture geometry that maximizes (or minimizes) ∆ ),- divided by the fracture area produced by fracture propagation, ∆, 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,  ., tensile strength,  ., and internal friction coefficient,  / .
Three of the scenarios use the same  .and  .with varying internal friction,  / , while the fourth uses higher  .and  .(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 ( / =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  !(Figure 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 Figure 4 (Figure S8).We refer to these different loading steps as the first and second stages of growth.For the models with  / =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 (Figure 4), but instead includes a diffuse distribution of fractures in a grid pattern (Figure S8, Figure 5e, f).This result indicates that fracture networks can 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  !(Figure 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 (Figure 5g, h).
For all the combinations of strength parameters, the fracture networks develop an increasing proportion of extensile/shear fractures, similar to the experiments (Figure 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.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,  .=5 MPa,  / =1.00 (Figure 5d, h), appears similar to the geometry of wing cracks that grow from inclined, precut fractures in experiments [44], and in numerical models [45].In contrast, lower intact rock strength promotes shear failure that links neighboring fractures in more planar geometries, and does not produce wing cracks (Figure 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 (Figure 3, Figure S4).After this 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 (Figure 3, Figure 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 sheardominated deformation under higher differential stress later in loading, and following the peak stress (Figure 1c-d).The observed evolution from shear to extensile deformation approaching failure agrees with machine learning analyses that use strain data from triaxial compression experiments and show that the local dilative strain provides more accurate predictions about the timing of macroscopic failure than the local shear strain [46], 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 [47,48].Microstructural observations acquired by serial thin-sectioning suggest that the first stress-induced fractures develop at preexisting intergranular and healed transgranular fractures in granite [22].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  !direction, consistent with our observations.
The evolution observed in the present analysis 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 [22].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 [4].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 of narrow columns separated by arrays of fractures aligned subparallel to the  !direction [12,49].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 [45].
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 [17].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 as similar to the friction between two surfaces because significant sliding may only occur on fractures following macroscopic failure [34].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 [50].However, experimental and numerical analyses suggest that the frictional sliding of fractures preceding macroscopic failure contributes to the internal friction coefficient [51], and thus the internal friction is some function of the strength of intact regions between fractures and frictional sliding on fractures [52].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 (Figure 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  ' at failure.
Numerical modeling provides insight into the experimental observations, and a physical interpretation of the evolution from shear to extensile-dominated deformation (Figure 5).
Constraining the geometry of first few fractures that fail during triaxial compression using energy optimization indicates that these fractures will not be aligned sub-parallel to the  !direction, but instead at 30-35° from the  !direction (Figure 4), consistent with the higher proportion of shear fractures observed early in loading in the experiments (Figure 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 [35].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 (Figure 5a-d, Figure 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 (Figure 5e, f, Figure 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 [53,54].Such delocalization is often attributed to the presence of heterogeneities [55].Several large earthquakes in southern and Baja California were preceded by the localization of low-magnitude seismicity that included phases of delocalization [56].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 (Figure 5i), consistent with the experimental observations (Figure 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.

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 [57] 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,  !, in steps of 0.5-5 MPa, with smaller steps closer to failure, until the rock failed with a stress drop (Table S1, Figure 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  !, we acquired an X-ray scan within 1.5 minutes 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 fracture development.We damaged the cores by heating them in an oven (Table S1).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 1600x1600x1600 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 TM , including the application of a non-local means filter [58].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 [59].This thresholding technique is not strongly influenced by noise [59] and produces segmented scans with porosity similar to values measured with sample imbibition [60].We then apply several processing techniques to the segmented scans in Avizo3D TM , 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  !direction.
Because fractures must have apertures greater than the spatial resolution of the tomograms (6.5 micrometer/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.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 (=0-17°) is similar or slightly larger than the average anisotropy of the shear-dominated fractures (=17-32°) throughout loading in each experiment (Figure S9).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 labelled 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 [61].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 [33], 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., 22].

Numerical model design
To identify the arrangement of fractures in the first stage of fracture growth, we use the software Fric2D [62].Fric2D uses continuum mechanics to calculate the stresses and strains within a homogeneous linear elastic material containing fractures and frictional faults in twodimensional plane strain [62].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  ' = 45 MPa (Figure 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 [27][28].
To simulate the growth of fractures, we use the software GROW [43], which calls Fric2D.
GROW accurately simulates the propagation path and linkage of two nearby fracture tips separated by a releasing step [63], the propagation of wing cracks from inclined fractures observed in uniaxial and biaxial compression experiments [46], the coalescence of hundreds of fractures observed in uniaxial compression experiments [64], and the development of thrust faults in accretionary prisms [65].
Recent work describes in detail the key algorithms of GROW [64,65].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.Table S2 and Table 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 [66].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 [3].For low porosity, nominally intact, crystalline rocks such as granite, the inherent shear strength ranges from 50-100 MPa, and the tensile strength ranges from about 5-20 MPa [3,67,68].Following the approach of Fattaruso et al. [64], 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  ., for constant  / and  ., 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 (Figure S10), consistent with the experiments.

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 [69].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  !direction [27][28].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 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  !under compression, where  is the angle of internal friction [23 p. 95-97].

Figure 1 .
Figure 1.a-b) An explanation of development of parallel and sub-parallel to the  !direction, to the volume of fractures aligned along the orientations expected for shear following the Mohr-Coulomb criterion, about 17-32° from the  !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,  & , of 5-20 MPa, fluid pressures,  ( , of zero to 10 MPa, and effective stresses, confining stress minus fluid pressure, of 5 MPa and 10 MPa.This suite of experiments provides an unparalleled 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.

Figure 2 .
Figure 2. Two-dimensional slices and three-dimensional volume renderings of fractures identified in X-ray tomography scans of experiments performed on intact (a) and damaged 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 sheardominated.The trends that we highlight for the combination in which extensile-dominated fractures have orientations,  ) =0-17°, and shear-dominated fractures have orientations,  * =17-32°, also occur for the other combinations (e.g., FigureS4-S5).

Figure 3 .
Figure 3. Evolution of the ratio of the volume of the extensile-dominated (=0-17°) and

Figure 4 .
Figure 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-

Figure 5 .
Figure 5. Fracture geometries propagated using the optimization of work with different