Catastrophic Failure: How and When? Insights From 4‐D In Situ X‐ray Microtomography

Catastrophic failure of brittle rocks is important in managing risk associated with system‐sized material failure. Such failure is caused by nucleation, growth, and coalescence of microcracks that spontaneously self‐organize along localized damage zones under compressive stress. Here we present X‐ray microtomography observations that elucidate the in situ micron‐scale processes, obtained from novel tri‐axial compression experiments conducted in a synchrotron. We examine the effect of microstructural heterogeneity in the starting material (Ailsa Craig microgranite; known for being virtually crack‐free) on crack network evolution and localization. To control for heterogeneity, we introduced a random nanoscale crack network into one sample by thermal stressing, leaving a second sample as‐received. By assessing the time‐dependent statistics of crack size and spatial distribution, we test the hypothesis that the degree of starting heterogeneity influences the order and predictability of the phase transition between intact and failed states. We show that this is indeed the case at the system‐scale. The initially more heterogeneous (heat‐treated) sample showed clear evidence for a second‐order transition: inverse power law acceleration in correlation length with a well‐defined singularity near failure and distinct changes in the scaling exponents. The more homogeneous (untreated) sample showed evidence for a first‐order transition: exponential increase in correlation length associated with distributed damage and unstable crack nucleation ahead of abrupt failure. In both cases, anisotropy in the initial porosity dictated the fault orientation, and system‐sized failure occurred when the correlation length approached the grain size. These results have significant implications for the predictability of catastrophic failure in different materials.


Introduction
Catastrophic failure of rocks in the brittle Earth is a critically important driving mechanism for phenomena such as landslides, volcanic eruptions, and earthquakes, including induced seismicity. Such failure often happens suddenly and with devastating consequences, occurring when structural damage, in the form of smaller faults and fractures, concentrates within localized zones. Damage localization leads to weakening and stress redistribution, eventually resulting in system-sized brittle failure along a distinct and emergent fault plane. Localized damage is pervasive at all scales throughout the brittle crust (Mitchell & Faulkner, 2012) and is therefore a fundamental control on catastrophic failure. Crack nucleation and growth, and crack coalescence within already localized zones, are relatively well-understood from microstructural and field observations of damaged rocks and from monitoring and locating earthquakes and acoustic emissions (elastic wave packets released during laboratory-scale microfracturing events). However, the process of localization remains elusive. Smaller cracks spontaneously self-organize along the incipient fault plane, often immediately before failure, but the precise mechanisms involved in this self-organization have yet to be determined. Open questions include (i) how do cracks, pores, and grain boundaries interact locally with the applied stress field to cause catastrophic failure to occur at a specific place, orientation, and time? (ii) Why can we detect precursors to catastrophic failure only in some cases? Fractures and faults have a self-similar structure; they are scale-invariant in their length and spatial distributions (Bonnet et al., 2001;Main, 1996;Main et al., 1990), and in the way their size relates to the energy released during rupture (Abercrombie, 1995;Goodfellow & Young, 2014). Remarkably, earthquakes and acoustic emissions (AE) are indistinguishable apart from the absolute source size, with scaling characteristics that are invariant over 15 orders of magnitude (Goodfellow & Young, 2014). This behavior is controlled by the local stress state and rock mass properties. Classically, brittle rock deformation and failure have been characterized by AE, with progressive cracking in heterogeneous materials under stress leading to systematic changes in the AE event rate and its frequency-magnitude distribution. Experiments have shown that pervasive microcracking accumulates in the sample subcritically, that is, without causing system-sized failure (Lei et al., 2000;Lockner et al., 1991), until the accumulating cracks self-organize along an asymmetric, localized damage zone. System-scale failure then occurs when nucleating microcracks have localized sufficiently for a runaway positive feedback of self-sustaining crack propagation and coalescence to take over (Main et al., 1993). In some cases, this self-organization becomes evident in the emergence of an inverse power law acceleration of event rate with a well-defined failure time. In others, system-sized failure of rock samples is commonly associated with the transition from an exponential increase to a sudden, rapid acceleration in the AE event rate close to peak stress (Moura et al., 2005;Sammonds et al., 1992;Vasseur et al., 2015). This transition occurs exactly when cracks begin to localize along the incipient fault plane (Lockner et al., 1991). At this crucial point, nucleated cracks grow by jumping geometrical and rheological barriers, so regions of stress concentration must already be correlated at the scale of the incipient fault network (Sammis & Sornette, 2002). The organized fracture network then propagates dynamically, with macroscopic failure occurring at a well-defined, finite time as the power law reaches its asymptote.
This behavior indicates a transition from pervasive but stable crack growth, controlled by the sample's microstructure, to an unstable regime of dynamic rupture along an organized fracture network, controlled by stress and fracture mechanics (Alava et al., 2008;Guéguen & Schubnel, 2003). The inverse power law transition can be described as a critical or second-order phase transition; a continuous transition from one state to another, during which the system becomes extremely susceptible to external factors. It is second-order if the first derivative of the free energy of the system (an entropy term) changes continuously as a function of the control parameter, for example, temperature (Stanley, 1971, Figure 2.6) or, in the case of a constant strain (or stress) rate rock deformation experiment, strain (or stress). This is associated with an inverse power law acceleration of the correlation length toward the critical point (Bruce & Wallace, 1989). At this point, strong correlations exist between all parts of the system (including at long-range), and many length scales become relevant (resulting in a self-similar structure and power law scaling), with events occurring at all relevant length scales (associated with broadband self-similarity of correlations). The transition to an inverse power law in the AE event rate, with its "finite time singularity" at failure, is also characteristic of a second-order or critical phase transition (Sammis & Sornette, 2002). If this occurs in the lead up to macroscopic failure, then the failure time can be forecast accurately (Vasseur et al., 2015(Vasseur et al., , 2017. Inverse power law acceleration to a well-defined failure time has also been seen in the evolving microstructure (microcrack porosity and the volume of the largest microfracture) of crystalline rocks undergoing brittle deformation (Renard et al., 2017(Renard et al., , 2018. system free energy. In turn, this depends on the applied stress, the length of the starting flaw and the specific surface energy of the material (Griffith, 1921(Griffith, , 1924. This results in a discontinuous or first-order transition between intact and failed states. In real materials that possess only a small amount of microstructural disorder, progressive subcritical cracking, that is, cracking which does not fulfill the conditions for sustained propagation, shows only an exponential increase in the event rate timeto-failure behavior (Vasseur et al., 2015(Vasseur et al., , 2017. Failure occurs suddenly and early; often much earlier than expected from an exponential asymptote (Vasseur et al., 2015(Vasseur et al., , 2017. This behavior is also characteristic of an abrupt first-order transition, with the exponential behavior reflecting local correlations (Sethna, 2006;Stanley, 1971).
Phase transitions are often characterized by the evolution of the correlation length and scaling exponents of the system in question (Stanley, 1971). The correlation length, ξ, is the distance over which the effects of a local perturbation of the system will spread (Thouless, 1989). Close to a critical point, the system can be viewed as made up of regions of size ξ in the critical state. In this case, ξ can be interpreted as the size of the regions of the nucleating phase, or the typical linear dimension of the largest piece of correlated spatial structure (Bruce & Wallace, 1989), which in our case is approximately the length of the largest fracture. As the two phases (intact and failed) compete to select the final equilibrium state, regions closer than ξ are correlated, and those further apart are uncorrelated. Approaching the critical point, the correlated (nucleating) regions become comparable to the system size. Thus, the maximum correlation length and associated parameters, such as the maximum AE magnitude, are restricted by the system size.
During a first-order transition, the correlation length ξ becomes macroscopically large but remains finite until the discontinuity at the sudden change of state (Stanley, 1971). In the case of a single Griffith crack, the correlation length is simply the length of the starting flaw, which suddenly becomes system-sized at failure as the flaw propagates instantaneously through the material (Figure 1, blue line). When real, structurally homogeneous materials, with a dilute population of material flaws, undergo progressive subcritical cracking (e.g., Vasseur et al., 2015Vasseur et al., , 2017, we expect the correlation length to increase exponentially but remain finite until it becomes system-sized at a sudden-onset discontinuity ( Figure 1, orange line). Conversely, during a continuous phase transition, the correlation length ξ becomes effectively infinite ( Figure 1, green line), growing as an inverse power law function of the control parameter (e.g., temperature), T, approaching the critical point, T c (Bruce & Wallace, 1989): This is the type of behavior associated with progressively more heterogeneous materials undergoing brittle deformation (Vasseur et al., 2015(Vasseur et al., , 2017. In summary, the degree of microstructural disorder of a material exerts a strong control on the type of phase transition from subcritical crack growth to dynamic rupture and consequently the predictability of the transition. In particular, experiments (Vasseur et al., 2017) and models (Kun et al., 2018) have shown that heterogeneity strongly influences the spatial distribution of microcracks at failure.
Here we investigate the impact of material heterogeneity on the nature of the phase transition between intact and failed states, and the associated predictability of failure, at the micron-scale. We show how the microcrack network evolves within a deforming crystalline rock with different amounts of disorder. Since pre-existing cracks are the most dominant factor of all heterogeneities that govern the fault nucleation process in laboratory rock samples (Lei et al., 2000), we deformed two samples of Ailsa Craig microgranite: one being an as-received control (nominally crack-free) and the other containing a pre-existing nanoscale crack network, induced by thermal stress, as a proxy for increased heterogeneity. Ailsa Craig samples, as received from the quarry, have no detectable cracks on thin sections under both optical and scanning electron microscopes (Meredith et al., 2005. They are an extreme end member of lowest crack density in natural rocks. Through the analysis of 4-D, in situ synchrotron X-ray microtomography (μCT) images of the two samples undergoing tri-axial deformation (see Cartwright-Taylor et al., 2020, for access to the data), we test the hypothesis that the transition to failure is abrupt and unpredictable (first-order) in the as-received sample (our initially crack-free end member) but is continuous and predictable (second-order) in the precracked sample. In situ observation of the deforming microstructure allows us to measure directly the relevant parameters such as the correlation length and the scaling exponents.
We find that increasing the microstructural disorder affects the geometry, size, and spatial distribution of the evolving microfractures. Using a combination of visual inspection of the μCT images, geometrical analysis of the evolving crack network, and techniques used in statistical seismology, we show that the microcrack network evolution varies significantly between the two samples. The degree of starting heterogeneity controls (i) the evolving spatial clustering and anisotropy of the microcracks and (ii) the order of the phase transition. The initially crack-free sample exhibits an exponential increase in damage that reflects local correlations, a finite correlation length, and no obvious precursors to failure. In contrast, the precracked sample exhibits emergent power law behavior, an inverse power law acceleration to infinite correlation length, and clear precursors to failure. However, though the parameters may be different, some crack growth characteristics and the power law scaling of the microcrack volume and intercrack length distributions appear independent of heterogeneity. Allowing for the fact that such microscopic failure characteristics may not be detectable above ambient noise in a field experiment, this may explain why measureable geophysical precursors to catastrophic failure events are detected only in some cases.

Ailsa Craig Microgranite and Thermal Stressing
We used Ailsa Craig microgranite (ACM) from the island of Ailsa Craig in the Firth of Clyde, Scotland. ACM is an extremely rare type of silica-unsaturated, alkali-rich microgranite, known as Blue Hone . As received from the quarry, no pre-existing microcracks are detectable either by optical or scanning electron microscopy (Meredith et al., 2005. Porphyritic in texture with a groundmass of mean grain size 250 μm, ACM contains sparsely distributed microphenocrysts (up to 1.5 mm) of alkali feldspar (Odling et al., 2007). Clint et al. (2001) found it to have extremely low porosity (<<1%) and permeability (1.5 × 10 −23 m 2 at 10 MPa effective pressure), indicating that the small amount of pre-existing pores are predominantly unconnected (e.g., see Figure 3 top left in Meredith et al., 2012). These properties make ACM ideal both for its main commercial use as the running surface of curling stones and for the purposes of this study. We chose ACM for two main reasons: (i) its small grain size (250 μm) and (ii) its virtually crack-free nature. The former is essential to ensure a statistically significant number of grains (>10 grains per diameter) in the small (3 mm diameter × 9 mm long) cylindrical samples and so to ensure that such small samples are representative of the rock as a whole. The latter is essential to allow comparison between two extreme end members: (i) an as-received control sample with the lowest possible (to our knowledge) pre-existing crack density and (ii) a second sample (from the same small block) containing a thermally induced crack network imprinted over the nominally crack-free microstructure, thus increasing its heterogeneity compared with the initially crack-free (untreated) sample.
To introduce a network of microcracks, one sample was heated slowly to 600°C prior to deformation. Thermal stressing is one of the key fracture-generating mechanisms in crustal rocks and is an effective method for introducing microfractures into rock samples. Heating ACM to elevated temperatures (>500°C) induces significant, permanent microcrack damage, evident from photomicrographs  and up to 50% and 30% reduction in P and S wave velocities, respectively (Clint et al., 2001). Scanning electron micrograph observations (Odling et al., 2007) show that heating ACM to 900°C causes the formation of a permanent microcrack network with average aperture of 0.3 μm formed by tensile failure, with each crack nucleating halfway along a previous one to generate fracture intersections of primarily T-shaped geometry. The thermally induced crack network is not discernible in our μCT data because this aperture is less than one tenth the length of one pixel (2.7 μm). Due to the partial volume effect, microcracks with an aperture smaller than half a pixel are not visible (e.g., Voorn et al., 2013).

Experimental Apparatus, Sample Assembly, and Loading Protocol
Synchrotron X-ray μCT, in combination with X-ray transparent pressure vessels (e.g., Butler et al., 2017;Fusseis et al., 2014;Renard et al., 2016), allows the microstructural evolution of deforming rock samples to be imaged directly, noninvasively, and in situ during an experiment. This provides a critical advantage over conventional deformation experiments, where the evolution of microscopic deformation cannot be inferred from post-test analysis of the microstructure because it is heavily overprinted by extensive damage caused during the macroscopic rupture process. Even in the case where conventional experiments are stopped immediately prior to macroscopic failure, overprinting occurs when the hydrostatic and differential stresses are released during extraction of the sample from the steel-walled pressure vessel, resulting in permanent damage and hysteresis. In situ X-ray μCT imaging overcomes both these issues, as well as providing detailed microstructural information about the temporal evolution of damage accumulation at a much higher temporal resolution. A single time-resolved experiment is equivalent to tens of conventional experiments with ex situ, post-experiment analysis, and has the virtue that the same sample is observed at each time step rather than a suite of samples, removing the issue of sample variability.
In our experiments, each sample of ACM underwent tri-axial deformation to failure. The experiments were conducted using a novel, miniature, lightweight (<1.3 kg) X-ray transparent tri-axial deformation apparatus 'Mjölnir', developed and tested at the University of Edinburgh. Mjölnir, named for the hammer of Thor, the Norse god of thunder, accommodates samples of 3 mm diameter and up to 10 mm in length and is designed to operate up to confining pressures of 50 MPa and axial stress in excess of 622 MPa (Butler et al., 2017). For this study, Mjölnir was installed on the μCT rotation stage at the PSICHE beamline at SOLEIL Synchrotron, Gif-sur-Yvette, France (Figures 2a and 2b). Two cylindrical samples of ACM, one heat-treated and one untreated, were cored using a diamond core drill and the ends ground flat and parallel to achieve 3 mm outside diameter and 9 mm length, compared to the typical grain size of 250 μm. Even though this sample diameter is very small (required to obtain high-resolution μCT images), the small grain size means that there are more than 10 grains per diameter, ensuring that such small samples are representative of the rock as a whole. The sample was assembled between the two pistons, jacketed with silicone tubing (3.18 mm internal diameter and 0.79 mm wall thickness), and protected from the confining fluid using twisted wire loops to seal the jacket against the piston (Butler et al., 2017). The pressure vessel was lowered over the sample assembly and fixed into place. A confining pressure of 15 MPa was then applied and maintained during the test. A hydrostatic starting pressure condition was achieved by simultaneously increasing the axial pressure to match the confining pressure. Delivery of the pressurizing fluid, deionized water, to the hydraulic actuator and pressure vessel was achieved using two Cetoni neMESYS ™ high pressure syringe pumps operated with QmixElements ™ software. Experiments were conducted at room temperature under nominally dry conditions. A reference μCT scan was acquired at zero differential stress to obtain the initial state of the sample prior to loading. The sample was then loaded to failure at a constant strain rate of 3 × 10 −5 s −1 in a step-wise manner, with steps of 20 MPa to start with, decreasing to 10 MPa from 70% of the failure strength and then 5 MPa once the sample started to yield (Figure 2c). At each step the stress was maintained, and a μCT volume acquired. To accommodate the full sample length at maximum resolution, three sequential scans were acquired at different positions along the length of the sample and then stacked. For each position the corresponding projections that comprised the full length of the sample were tessellated and merged to create a single projection used for reconstruction of the whole sample in one μCT volume. Each full set of scans was acquired in approximately 10 min. For each sample, 15 sets were acquired during loading with an additional set acquired after the main failure. For the heat-treated sample, this included one set at peak differential stress of 185 MPa. This μCT volume contained the incipient fault at the critical point of failure, and the sample failed immediately upon continuation of the loading procedure. The untreated sample reached a peak stress of 182 MPa but failed before it could be scanned at this stress. The last pre-failure scan was at 177 MPa (97% of the critical failure stress, σ c ).
The differential stress is σ = σ 1 − σ 3 , where σ 1 is the axial stress (the product of the measured ram pressure and the difference in area between the ram and the sample cross-section) and σ 3 is the radially symmetric confining pressure. Axial sample strain was calculated as ϵ = δL/L 0 , where δL is the change in length of the sample between the starting μCT volume and the volume of interest and L 0 is the initial length of the sample. It was obtained directly from the reconstructed μCT volumes by measuring the length change of the rock core between two fixed locations in each volume.

X-ray Imaging and Image Data Preprocessing
X-ray μCT data were acquired using an sCMOS (scientific Complementary-Metal-Oxide Semiconductor) Hamamatsu ORCA Flash4.0 camera, with a Rodagon 50 lens, giving about 2.5× magnification (effective pixel size 2.7 μm) and a 250 μm thick LuAG:Ce scintillator. The white beam with an average detected energy of about 66 keV was filtered with 1 mm aluminum and 0.5 mm tungsten. During each scan, 1,200 projections were acquired over 180°, with an exposure time per projection of 15-19 ms depending on the progressive darkening of the objective lens. A mix of absorption and phase contrast data was acquired, with a sample to detector distance of 125 mm.
Each μCT volume was reconstructed by filtered back projection. Reconstructions were performed at the PSICHE beamline, using both X-ray absorption and phase contrast modes as implemented in the PyHST2 software (Mirone et al., 2014), and yielded 3-D volumes of 1,700 × 1,700 × 4,102 equidimensional voxels, with a voxel edge length of 2.7 μm. These volumes were then processed to extract the fracture network from the reconstructed images. To deal efficiently with the huge size of each 3-D volume (approximately 40 GB) and speed up the subsequent processing, we selected a subvolume of interest-the region in the failed samples where the majority of damage had accumulated ( Figure 2d, Table 1). Using the Avizo ™ software package, this subvolume was extracted from each of the full 3-D volumes and de-noised with an anisotropic diffusion filter (stop value 0.4 over four iterations), which was chosen to emphasize the microcrack features as it preserves strong edges and enhances edge contrast. It was then down-sampled to 16-bit with a 32-bit threshold range of −0.3 to 0.8, yielding individual data sets of manageable size (approximately 3 GB).

Segmentation
Each data set of 3-D gray-scale images was then segmented to separate from the rest of the rock matrix the pre-existing pores and the evolving deformation-induced cracks in a binary fashion. Herein we use the term "porosity" to include all the segmented void space in the sample, whether pre-existing (and therefore associated with the igneous history of the rock) or deformation-induced. We use the term "void" to describe an individual segmented object.
Although easily distinguishable by the human eye, narrow planar features such as fractures are difficult to extract automatically from large 3-D image data sets. This is due to the range of gray-scale values accommodated by fractures of different apertures and the increasing similarity of these gray values to the surrounding rock matrix as the aperture decreases. The main reason for this is the partial volume effect, whereby voxels containing both air and rock matrix appear brighter than voxels containing air alone. Fracture surface roughness and narrow apertures contribute to this effect. We used a multiscale Hessian fracture filter (MSHFF) technique to meet these challenges while still using an automated approach and segment the microcracks from the image data. This technique, developed and described in detail by Voorn et al. (2013), uses the Hessian matrix (second-order partial derivative of the input image data) to represent the local curvature of intensity variation around each voxel in a 3-D volume (e.g., Descoteaux et al., 2005). Attributes of this local curvature can be used to distinguish planar features in the data set (Text S1a in the supporting information). The analysis is conducted over a range of observed crack apertures, which are combined to produce the final multiscale output: narrow fractures of varying apertures detected within the 3-D image data. The analysis was carried out using the macros for FIJI (Schindelin et al., 2012) published by Voorn et al. (2013), utilizing the FeatureJ plugin (Meijering, 2010) to calculate the Hessian matrices, with input parameters given in (Table 2).
The binary segmented volumes were labeled, and an intensity radius threshold of 15 to 100 was applied followed by a small spot filter to remove round segmented objects with a radius of <50 pixels (134 μm) for visual clarity of the fault plane ( Figure 3a). The local aperture at each voxel of the segmented voids was computed from the diameter of the largest ball containing the voxel and entirely inscribed within the object (Hildebrand & Rüegsegger, 1997). Even with the segmentation method described, there is still significant undersampling of the void population, particularly at the narrower end of the aperture range ( Figure 3b). Further work in this area is required and would benefit from machine learning approaches (Andrew, 2018). Note. Definitions given in Voorn et al. (2013). We present the data according to a co-ordinate system (x, y, z) where z is the vertical axis, which is parallel to the loading direction and corresponds to the direction of axial stress (σ 1 ). The other two (x and y) are the horizontal axes, which are perpendicular to the loading direction and correspond to the confining stress (σ 2 = σ 3 ) with their directions arbitrarily assigned but consistent between the two experiments. Void orientations are given in terms of their dip ϕ (deviation from horizontal) and strike θ (deviation from y).

Analysis of the Segmented Microcrack Network 2.5.1. Porosity and the Number of Fractures
Initially of interest is the evolution of porosity, φ (including both the pre-existing pores and the developing microcracks) and the total number of voids, N, in each 3-D subvolume as deformation progressed. Due to the irregular shape of the segmented objects, and to take into account crack coalescence, one void was defined as all objects that were connected at least by the apex of a corner. To obtain N, each individual void, i, was assigned a label using the Label Analysis module in Avizo ™ . Porosity was obtained from the total void volume divided by the analyzed subvolume: φ = (∑V i )/(πR 2 l), where V i is volume of each crack, R is the sample radius, and l is the length of the analyzed subvolume. To determine the most likely empirical relationship with strain for both φ and N, we found the parameters for several possible models (quadratic, exponential, and simple power law) using nonlinear least squares regression and then used the corrected Akaike information criterion, AICc (Hurvich & Tsai, 1989), to test these competing models objectively (see Text S2 for a full description of the calculations).

3-D Microcrack Orientations and Geometries
These were obtained from the binary segmented data using an object-based approach to determine the best-fitting ellipsoid around each segmented void (pore or microcrack). Each ellipsoid was calculated independently from the crack's 3-D moments of inertia (Text S1b), also using the Label Analysis module in Avizo ™ . First-order moments, M, define the void's center of mass (centroid). Second-order moments, I, define the inertia (or covariance) matrix, with eigenvectors representing the ellipsoid axes orientations.
Major, minor, and medium ellipsoid radii, r, were computed as r¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 5 eigenvalue j j p from each corresponding eigenvalue of the inertia matrix (Ollion et al., 2013).

Crack Network Scaling Exponents and Correlation Length
To find evidence for the type of phase transition undergone by each sample, we obtained the following indicators of critical point behavior: the correlation length, ξ (linear dimension of the largest void), the void size exponent, β, and the void separation exponent, or correlation dimension, D. These time-dependent parameters are equivalent to those commonly used to quantify the evolution of seismicity from acoustic monitoring (Aki, 1965;Bufe & Varnes, 1993;Kagan, 2002;Ouillon & Sornette, 2000;Sammis & Sornette, 2002;Sornette & Sammis, 1995;Sykes & Jaumé, 1990;Turcotte, 1997;Tyupkin & Giovambattista, 2005;Zöller et al., 2001). In rock deformation studies (e.g., Lei & Satoh, 2007;Moura et al., 2005), these parameters have been similarly inferred from the amplitudes and locations of acoustic emissions (AE). In particular, the intercrack distance inferred from AE is a key parameter that controls the failure time and hence the accuracy of failure time forecasts (Vasseur et al., 2017).
In this study, we obtained ξ, β, and D directly from the evolving population of microcracks in the μCT data, rather than indirectly from AE. We used void volume as a metric for void size and estimated inter-void distances (void separation) from the distribution of points defined by ellipse centroids. We obtained void volumes, V i , and centroids (Text S1b) directly from the Label Analysis module in Avizo ™ and then computed Euclidean lengths, L i , between centroids.
We obtained maximum likelihood estimates of the void size exponent, β, from the frequency-volume distribution in each μCT subvolume. We tested three different but related models often used to describe the seismic moment distribution in seismicity (Kagan, 2002; full details of this procedure are given in Text S3): i) GR: the Pareto distribution (a pure power law, equivalent to the Gutenberg-Richter distribution) with cumulative complementary (survivor) function where V i is the volume of each individual void and V max is the upper bound (maximum) void volume in the distribution. ii) TRP: the truncated Pareto distribution (similar to the GR but showing a power law taper in the tail where V c is the tapering corner volume of the distribution.
iii) TAP: the tapered Pareto distribution (equivalent to the modified Gutenberg-Richter relation which shows an exponential taper in the tail toward We defined a correlation length, ξ¼ ffiffiffiffiffiffiffiffiffiffi V max 3 p for the pure power law model or ffiffiffiffiffi ffi V c 3 p for the truncated models. The completeness volume, V t , is the smallest void volume at which 100% of voids in a space-time volume are detected (Mignan & Woessner, 2012;Rydelek & Sacks, 1989;Woessner & Wiemer, 2005) and is equivalent to the threshold of completeness in seismicity data. We obtained V t from the maximum curvature method (Roberts et al., 2015), that is, from the peak of the incremental frequency-volume distribution. This method is appropriate for the sharp-peaked distributions seen in our data ( Figure S1). In both samples, the number of voids with V ≥ V t exceeded 200, which is the minimum catalogue size required for reliable estimation of β (Roberts et al., 2015). We used a modified Bayesian information criterion (BIC) (Bell et al., 2013a;Leonard & Hsu, 1999) to find the most appropriate model (see Text S3 for full details of the calculation and Figure S2 for the full results), and thus obtained the most likely values of β and ξ for each distribution. The BIC is more appropriate for distributions of large sample populations investigated here than the AIC, and the results can be compared more directly with previous work (Bell et al., 2013a).
To distinguish the type of phase transition, that is, whether or not ξ followed an inverse power law acceleration to failure, we fit the data by nonlinear least squares regression to an inverse power law model of the form: y = k(x p − x) −p , where x p is the predicted value of the control parameter x at failure, k is a scaling factor, and p parameterizes the rate acceleration of y, all determined by nonlinear regression. The point of failure, x p , is defined by a mathematical singularity as y → ∞. This is directly analogous to the approach to a critical point in a second-order phase transition for the correlation length (Equation 1). It is also equivalent to material failure forecasting approaches based on the Time-Reversed Omori Law for aftershock decay (e.g., Bell et al., 2013b;Cornelius & Voight, 1994;Kilburn, 2003;Kilburn & Voight, 1998;Smith et al., 2009;Utsu et al., 1995;Vasseur et al., 2015Vasseur et al., , 2017Voight, 1988). We used stress as the control parameter instead of time because the stepped loading procedure that we conducted precludes realistic temporal rate estimates. Importantly, this model makes no a priori assumptions about any of the parameters. The predicted failure stress, σ p , is what would be available in real time, rather than the observed failure stress, σ c . By estimating σ p independently, we can quantify any systematic error in its estimation by comparing it to the observed failure stress and hence quantify any bias in a potential forecasting scenario. We used a trust region algorithm to minimize the residual sum of squares between the observed data and the model (see Text S5 for details). We also tested an exponential model y = h exp(qx) as an alternate hypothesis. We cannot use a simple criterion such as r 2 alone to determine the relative goodness of fit because the competing hypotheses have different degrees of freedom, so we used the AICc (see Text S2 for details of the calculations). It is based on the residual sum of squares and is considered more robust than r 2 alone because it takes into account the number of parameters in the model, penalizes models with more parameters, and can be used to determine the relative likelihood of the models given the data.
We obtained the two-point correlation dimension, D, from the relation P( where P(L i ) is the incremental probability distribution of the inter-void lengths, L i (see Text S4 for more detail). The exponent, D inc = D − 1, of P(L i ) in the identified power law region, 30 < L i < 1,350 μm, was obtained from a linear regression in log-log space ( Figure S3). D is then D inc + 1. If D < 3, voids are clustered spatially, but as D → 3 they become volume-filling and therefore less clustered (Robertson et al., 1995).

Results: Microfracture Network Evolution
Here we present our 4-D observations of the evolving segmented fracture network in each sample, together with the distributions of microcrack orientations. This is followed by results from our quantitative analysis of the fracture network. We first show the influence of heterogeneity on the evolution with strain of (i) stress, porosity, and the number of voids and (ii) void geometry, which reveals how the initial, small anisotropy of the crack network increased in the lead-up to failure. Following this, we test our hypothesis regarding the type of phase transition undergone by each sample by showing the transition to failure of the correlation length as a function of stress. Finally we present the evolution with strain of the correlation length and the scaling exponents of void size and separation. For the purposes of testing our hypothesis and for clarity of presentation we have analyzed data up to failure but not beyond.

4-D Observations of Microfractures and Their Orientations
The spatial evolution of microcracks differs significantly between the two samples (Figures 4 and 5).
Although the untreated sample appeared to fail along a localized shear fault (Figure 4p), pre-failure damage accumulated in a spalling pattern of radial damage zones subparallel to σ 1 (Figures 4f-4o). Conversely, the heat-treated sample failed along a localized shear fault, inclined at 30°to σ 1 , along which pre-failure damage had already accumulated (Figures 5k-5p). In both samples, the localized damage zones consisted mainly of shear and axial microfractures oriented between 5°and 30°to σ 1 (Figures 6a and 6b) with bridging ligaments. Local fracture aperture increased as fractures propagated within the localized zones of damage.
In the untreated sample, damage localization (established visually from the segmented CT volumes in Figure 4) first occurred along four narrow zones simultaneously (orange ellipses in Figure 4f) at 1.37%  Table 3. Color coding shows the local fracture aperture at each voxel. See Data Set S1 for the zip file of high-resolution images.
strain and 64% of the failure stress, σ c , (Table 3). Damage was concentrated along these zones until further radial zones developed on the other side of the sample (pink ellipses in Figure 4j) at 85% σ c . It is not clear exactly how the eventual fault developed because the sample failed before the next scheduled image time.
The post-failure μCT subvolume (Figure 4p) indicates that slip occurred along a plane not previously localized and located slightly above where the radial damage had accumulated, which is consistent with the sudden, abrupt nature of the failure inferred from the stress-strain data (section 3.3.1). However, close to peak stress (97% σ c ; Figure 4o) cracks had begun to localize along a shear zone that was above and formed from the tip of some of the radial damage zones and was conjugate to the eventual fault plane.
In the heat-treated sample, damage localization (established visually from the segmented CT volumes in Figure 5) first occurred along a shear zone conjugate to the eventual fault plane (orange ellipse Figure 5h) at 1.24% strain and 72% σ c (Table 4). Damage progressively concentrated along this plane until  Table 4. Color coding shows the local fracture aperture at each voxel. See Data Set S2 for the zip file of high-resolution images.
localization around the incipient fault plane became apparent (pink ellipse in Figure 5k) at 90% σ c and 1.38% strain; the same amount of strain as the initial localization in the untreated sample. Fracture nucleation and propagation within the initial damage zone then stopped, continuing instead along this more favorably oriented zone until failure. This flip in orientation between two optimally oriented, conjugate, shear planes has previously been seen by Lennartz-Sassinek et al. (2014). Here it coincides with reduced sample stiffness and strain hardening inferred from the stress-strain data (section 3.3.1).
Microcrack dips, ϕ, became progressively more vertical with increasing stress in both samples (Figures 6a  and 6b), indicating the preferential nucleation of tensile cracks. These cracks formed en-echelon (Kranz, 1979;Olson & Pollard, 1991;Reches & Lockner, 1994;Tapponnier & Brace, 1976) and wing-crack (Ashby & Hallam, 1986;Ashby & Sammis, 1990;Fairhurst & Cook, 1966;Horii & Nemat-Nasser, 1985Nemat-Nasser & Horii, 1982;Nemat-Nasser & Obata, 1988) arrays (Figure 7), concentrated in the heat-treated case at the tip of the propagating fault zone. All radial damage zones in the untreated sample grew in this manner immediately after their initial localization (Figure 4f onward). In contrast, this process occurred only in the heat-treated sample during localization around the eventual fault plane (Figure 5l onward), not during the initial localization around the unfavorable conjugate. En-echelon and wing-crack arrays formed at 1.45 ± 0.01% strain in both samples (Figures 4g and 5l). At this point the untreated sample was only at 70% σ c , compared to 90% for the heat-treated sample.
One advantage of the 3-D sampling enabled by μCT imaging is that we can test the null hypothesis that the initial sample porosity is isotropic. The optimal strikes of the segmented voids, θ, show a predominant orientation in the initial porosity in both samples (Figures 6c and 6d). This starting anisotropy was more pronounced in the heat-treated sample than in the untreated sample (33.0 ± 15.1% compared with 14.3 ± 11.8%; see Table S1). Overall, anisotropy in the void strike increased steadily throughout  Table S1). The strike of the eventual fault closely followed this pre-existing anisotropy in both samples but to a much greater degree in the heat-treated sample. In the untreated sample, although the distribution peaks and troughs broaden as the radial zones localized, the strike of the post-failure fault was oriented within 30°of the initial preferred strike orientation (Figure 6c). In the heat-treated sample, the strike of the emerging fault plane tracked the orientation of the initial crack porosity anisotropy almost exactly (Figure 6d), while the distribution of peaks and troughs remained stationary and became more defined, as the damage zone localized.
Stereonet depictions (Figure 8) of the void orientations (poles to planes) projected along the axial direction (σ 1 ) confirm these observations, showing a predominant strike parallel to the pre-existing porosity in both samples, followed by the development of mainly vertical cracks at localization, in line with our visual examination. These cracks initially cluster along the pre-existing strike in both samples but become increasingly distributed in the untreated sample (blue stereonets) during yield and approaching failure, with failure occurring along a fault offset by around 30°to the pre-existing strike. Conversely, in the heat-treated sample (orange stereonets), these vertical cracks cluster increasingly along the pre-existing strike throughout deformation. Closure of some shallow dipping voids is seen in both samples.
In the time-lapse video of the untreated sample in the x, y projection (Video S1), the first axial fracture initiated at a spot on the sample edge close to the bottom of the sample, below and on the outside of a region of concentrated porosity. Further localization occurred simultaneously along vertical zones distributed radially around the sample, which appear to have grown preferentially into the sample before propagating vertically as crack segments linked up. Approaching peak stress, an array of microcracks with varying orientations formed around the region of concentrated initial porosity, bridging three radial zones in a curved damage zone. The same process occurred again at 97% σ c , bridging four radial zones adjacent to the previous three but in a conjugate orientation. This bridging fault propagated up the sample (Video S2; y, z projection), at a different strike to the post-failure fault.
Time-lapse video of deformation in the heat-treated sample in the x,y projection (Video S3) shows that localization initiated within the sample, not at the boundary, on the site of a pre-existing void, precisely as anticipated by Griffith (1921Griffith ( , 1924. Subsequent microcracks that localized along the damage zone nucleated between the initial site and the sample edge. The emerging fault plane initiated at the sample boundary and grew horizontally into the sample, as subsequent microcracks localized along it, before propagating down the sample parallel to the z, x plane (Video S4; y, z projection). Simultaneously, microcracks localized on the opposite side of the sample along the same strike as the initial, arrested damage zone. As deformation continued and the sample reached peak stress, microcracks nucleated in the center of the sample. These cracks joined the optimally oriented damage zone to the conjugate damage zone on the other side of the sample, resulting in a curved shear zone, consisting of arrays of microcrack segments linked by bridges of intact rock, along which the sample failed. Note. Letters A-P refer to the image volumes shown in Figure 4. Localization first appeared in scan F along several vertical zones distributed radially around the sample. Additional zones localized in scan J, but scan P shows that failure occurred along an unrelated shear fault. Note. Letters A-P refer to the image volumes shown in Figure 5. Localization was first seen in scan H along a shear zone that subsequently arrested when microcracks localized along a more favorably oriented damage zone in scan K, along which the sample eventually failed.

10.1029/2020JB019642
Journal of Geophysical Research: Solid Earth

Evolution of Specimen and Microcrack Characteristics With Strain 3.2.1. Stress, Porosity, and the Number of Voids
Both samples had a similar stiffness early on (Figure 9a), but the heat-treated sample became stiffer at 0.84% strain when the number of voids decreased slightly (Figure 9c). This is consistent with compaction of the compliant thermal microcracks. The onset of localization, as determined visually from the µCT volumes, is evident in both samples as a yield point in the stress-strain curve; at 1.37% strain and 1.24% in the untreated ( Figure 4f) and heat-treated (Figure 5h) samples respectively. Further yielding occurred once the damage zone propagated at 95% σ c in the heat-treated sample (Figure 5m), but only from 97% σ c in the untreated sample (Figure 4o).
The heat-treated sample had lower pre-existing porosity than the untreated sample (φ 0 HT = 0.62φ 0 UT ) and fewer but slightly larger voids (N 0 HT = 0.54N 0 UT ), with half the number of voids accounting for two thirds of the porosity seen in the untreated sample (Figures 9b and 9c). However, this observation only accounts for voids visible above the detection threshold of the segmentation algorithm (a void Figure 7. Slice through the sample in the y,z projection showing in (a) the untreated sample close to peak stress ( Figure 4o) and (b) the heat-treated sample at peak stress ( Figure 5o) the formation of tensile (c) en-echelon (Olson & Pollard, 1991;Thomas & Pollard, 1993) and (d) wing-cracks (Fairhurst & Cook, 1966; Figure 15 in Damjanac & Fairhurst, 2010).

10.1029/2020JB019642
Journal of Geophysical Research: Solid Earth volume of 3,000-4,000 μm 3 ; see section 3.3.1) and does not include unresolved nanoscale thermally induced cracks. The observed differences may be accounted for by natural sample variation within these very small samples and/or some void closure from thermal expansion during the heat-treatment.
Both samples showed a tenfold increase in porosity, φ, over the duration of their respective deformation experiment (Figure 9b), but only a twofold increase in the total number of voids, N, in the heat-treated sample, compared with a nearly threefold increase in the untreated sample ( Figure 9c). This indicates that crack nucleation was more dominant in the untreated sample, compared with crack growth in the heat-treated sample. The untreated sample showed no evidence of compaction in the early stages of deformation and the onset of localization (Figures 4e-4f) is evident as a large jump in N of 600 voids at 1.37% strain and a corresponding threefold increase in φ (Figure 9). Conversely, in the heat-treated sample a small decrease in N of approximately 50 voids provides evidence for some early compaction due to void closure, although this equates to only a tiny proportion (0.005%) of φ 0 . This was associated with the closure of some optimally oriented (shallow dipping) voids prior to localization (Figure 8; orange stereonets). The onset of localization is evident as a minimum in both φ and N at 1.24% strain (Figure 5h), and both variables exceeded their initial values when the optimally oriented damage zone localized (Figure 5k). Once localization initiated, both samples showed an overall acceleration toward failure in both φ and N. However, in the untreated sample there were two occasions where the acceleration was temporarily arrested. The first of these corresponded to the propagation of new localized zones (Figure 4j), while the second corresponded to the change in orientation of the bridging zone (described in section 3.1). The heat-treated sample showed a slight slow-down in acceleration that corresponded to the nucleation of new microcracks between the two ends of the eventual fault (described in section 3.1), followed by a final acceleration immediately before dynamic rupture.
In both samples the evolution of both φ and N with strain is best described with simple power law models (Figures 9b and 9c); that is, they have the lowest AICc (Tables S2 and S3). The exponent is 3.1 for both Figure 8. Stereonet projection of void ellipsoid orientations (poles to planes) projected down the loading axis (i.e., in the x,y plane) for both samples. Letters refer to the stress-strain steps for the untreated (Table 3 and Figure 4) and the heat-treated (Table 4 and Figure 5) samples. Clustering contours were calculated from uniform kernel density estimation to show significant departures from a uniform distribution (Kamb, 1959). Kernel radius, r¼3= ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi π 9 þ n ð Þ p , where n is the number of data points. Contour intervals are given as a multiple of the standard deviation to emphasize the statistical significance of the number of points falling into each kernel (Haneberg, 2004). In the untreated sample, strikes were dominated by the pre-existing porosity until scan F, when microcrack localization initiated along the steeply dipping, radially distributed zones (orange ellipses in Figure 4). The radial pattern became increasingly symmetrical around the sample throughout the rest of the experiment as these zones propagated and as microcracks localized along new zones (pink ellipses in Figure 4). In the heat-treated sample, strikes were dominated by the pre-existing porosity throughout deformation, although the initial localized damage zone observed at scan H (orange ellipse in Figure 5) was conjugate to the eventual fault zone seen at scan J (pink ellipse in Figure 5).

Journal of Geophysical Research: Solid Earth
variables in the untreated sample, compared with 8.8 and 7.7 for φ and N respectively in the heat-treated sample, showing an acceleration toward failure that was almost three times faster in the heat-treated sample than the untreated one. These exponents also show that the acceleration in N accounted for all of the acceleration in φ in the untreated sample, confirming that crack nucleation was the dominant damage mechanism throughout deformation, whereas in the heat-treated sample, the acceleration in N did not completely account for all of the acceleration in φ, confirming that crack growth played an increasingly important role closer to failure.

Microcrack Geometry
To establish empirically how the microcrack geometry evolves with increasing deformation, we present the variation with strain of the mean value of the major, minor, and medium ellipsoid radii from the population of voids in each μCT subvolume (Figures 10a and 10b). We also show the mean ellipsoid eigenvalue ratios, used to infer the evolution of void aspect ratio (Figures 10c and 10d). We present two ratios: (i) the smallest to the medium eigenvalue of the covariance matrix (section 2.5.2 and Text S1b), where flatter objects have smaller values and (ii) the medium to the largest eigenvalue, where more elongated objects have smaller values.
Corresponding mean void radii, r, were about the same size in both samples. In the untreated sample (Figure 10a) r major and r medium (blue and orange circles, respectively) began to increase at the onset of localization (Figure 4f). In the heat-treated sample (Figure 10b) they began to increase as microcracks localized along the optimally oriented damage zone (Figure 5k), after the onset of localization. In both samples, r major were oriented approximately parallel to the strike of the eventual fault plane (Figures 6c and 6d), with their mean values, r major , increasing more quickly than r medium , showing that microcracks grew twice as fast along strike (i.e., perpendicular to σ 1 ) than down dip, becoming more elongate as deformation progressed (Figure 10c). Voids in the heat-treated sample were marginally flatter than those in the untreated sample, with voids in both samples becoming flatter as failure approached (Figure 10d). This implies that the scaling of crack growth, while scale-invariant in length, may be self-affine (variable aspect ratio) rather than self-similar (constant aspect ratio). The down dip extent (r medium ) increased from 2.5 to 4 times the crack aperture (r minor ; yellow circles). In the heat-treated sample, this was due to continued crack growth down dip relative to a constant crack aperture (Figure 10b). However, in the untreated sample, crack growth down dip stopped altogether close to failure and a small decrease in aperture accounted for the voids becoming flatter (Figure 10a). Growth along strike also stopped (within error), and the continued increase in the number of cracks at this stage ( Figure 9c) confirms that nucleation of new cracks accounted for almost all the porosity generation close to failure in this sample.
In summary, the average growth pattern of individual microcracks is independent of heterogeneity in the early stages of localization. The behavior changes close to failure, with continued crack growth in the heat-treated sample accounting for the faster acceleration in porosity than void number, while crack growth in the untreated sample was effectively halted in favor of crack nucleation, which accounted for the entire acceleration in porosity and void number as failure approached.

Evidence for Phase Transition Style
To establish the type of phase transition undergone by each sample, we present the evolution to failure of the correlation length and the scaling relations as a function of both differential stress σ and axial sample strain Figure 9. Evolution of (a) differential stress, σ, (b) porosity, φ, and (c) the number of segmented voids, N, with increasing axial strain, ϵ for the untreated sample (blue circles) and the heat-treated sample (orange circles). Dash-dot lines show the strain at which each sample failed, while dashed lines show the onset of damage localization as seen the μCT volumes (Figures 4 and 5). Solid lines in (a) show the region of data used to calculate the Young's moduli (11.3 and 9.5 GPa for the untreated and heat-treated samples respectively). Solid lines in (b) and (c) show the preferred simple power law models with 95% confidence intervals shown as dashed lines; see text for exponents. ϵ. Renard et al. (2018) argue that stress is a stronger control variable than strain, but strain is usually the only directly observable control parameter in real Earth applications. We first present the scaling relationships for void volume and inter-void length and then show how the correlation length, ξ (linear dimension of the largest void), evolves as a function of stress. We then analyze the evolution of ξ, β, and D (the void volume and inter-void length exponents respectively, defined in section 2.5.3) as a function of strain.

Microcrack Volume and Intercrack Length Distributions
Both samples show an approximately power law complementary probability distribution in void volume, V i (Figures 11a and 11b), with the proportion of larger voids increasing systematically with respect to strain and stress. Both samples also show an approximately power law distribution in their inter-void lengths, L i (Figures 11c and 11d), within a finite range, identified as 30 < L i < 1,350 μm (close to half the sample diameter), with little apparent difference in the shape of the distributions as stress and strain increase. We can therefore define power law scaling exponents β from the frequency-volume distributions and the correlation dimension D from the inter-void length distributions.
Values of completeness volume, V t (defined in section 2.5.3), ranged from 3,000 to 4,000 μm 3 , roughly equivalent to a void aperture of 14-16 μm. This is much larger than the theoretical detection threshold of half the pixel size (1.3 μm) consistent with undersampling of very narrow cracks during segmentation. Void volumes in the untreated sample are best described (i.e., have the lowest BIC) by the truncated Pareto distribution (TRP) at the three earliest steps of deformation and then by the characteristic Pareto distribution (GR), with the transition between the two models occurring at 43% σ c , two stages before the onset of localization ( Figure S2). In contrast, void volumes in the heat-treated sample are best described by the GR distribution throughout the experiment.

An Inverse Power Law Acceleration to Failure?
Parameters for an inverse power law acceleration to failure for the normalized correlation length, ξ/l, were obtained for both samples (Figure 12), using data observed in segmented μCT volumes between Figures 4f-4o and 5h-5o and the method described in section 2.5.3. While an inverse power law acceleration is commonly only distinguishable from an exponential acceleration within 10% of the singularity (Bell et al., 2013b), restricting the data to this region (stages L-O in the untreated sample and K-O in the heat-treated sample) would have left very few data points for the analysis. In the untreated sample the exponential and inverse power law models are indistinguishable over the data range, and the predicted failure stress, σ p , is far from the observed failure stress, σ c . The likelihood that the inverse power law model fits the data as well as the exponential model is just 3% (Table S5). Thus, it is impossible to define an accurate failure point in this sample. The sample failed abruptly, long before the predicted singularity, after an exponential acceleration in the correlation length. Conversely, in the heat-treated sample, σ p is accurate to within 3% of σ c , while the asymptote of the exponential model is further from σ c .The likelihood of the exponential model fitting the data as well as the inverse power law is 72% (Table S5), which, although relatively high, is not significant (>95%). These differences are diagnostic of a first (abrupt) and second (continuous) order phase transition respectively (Figure 1), validating our hypothesis that, within the temporal resolution of our experiments, the transition to failure is first-order in the untreated sample and second-order in the heat-treated sample.

Evolution of Crack Population Metrics With Respect to Strain
Both samples show a systematic increase in the correlation length ξ toward failure as a function of strain, with failure occurring when ξ increased beyond 200 μm (Figure 13a). This limit marks the longest crack supported by the sample volume without a runaway instability developing and falls just short of the mean grain size (i.e., the length between grain boundaries) of the groundmass (250 μm; see section 2.1). This implies that  Figure S1) using the models of Kagan (2002) and fitted the incremental L i data using linear regression in log-log space (see Figure S3) to find D, after Turcotte (1997). Blue distributions overlaying the rest are post-failure.

Journal of Geophysical Research: Solid Earth
the sample breaks when whole grains break. The nature of the increase in ξ(ϵ) in the untreated sample (blue) was exponential (i.e., had the lowest AICc; see Table S4) up to a finite ξ that fluctuated around 200 μm before failure. Conversely, in the heat-treated sample (orange) ξ(ϵ) preferred a simple power law acceleration ( Figure S4 and Table S4) to failure, failing when ξ > 200 μm. The exponent (6.9) is the same (within error) as the exponents for the evolution of φ and N with strain (section 3.3.1; Tables S2-S4), independently confirming that in this sample crack growth played an increasingly important role closer to failure. The power law acceleration emerged only once fractures began to localize along the optimally oriented damage zone (at 90% σ c and ξ > ξ 0 ), representing a strong self-organization in the crack network over all distances to concentrate on the damage zone that controls the eventual fault plane.
The two samples had different initial exponents β for their volume distributions ( Figure 13b). In the untreated sample (blue) β rose sharply to a peak at the transition from the TRP to the GR model. This shows that, at this point, the largest cracks in the taper were growing or opening, while simultaneously many more small cracks were becoming active above the segmentation detection threshold. The number of voids and the porosity were constant in this phase, implying other voids were simultaneously closing in compaction (Figures 9b and 9c). This trade-off is consistent with independent observations from acoustic emissions (Graham et al., 2010) and models (e.g., Brantut et al., 2012Brantut et al., , 2014 of the competition between compaction and dilatancy during the quasi-elastic region of the stress-strain curve (Figure 9a). Beyond this peak, β decreased smoothly to the first of two local minima once the additional radial zones had localized ( Figure 4k). This indicates instability in the sample-related crack nucleation (Figures 9c and 10a) and might be considered a precursor to failure, albeit without evidence of quasi-static damage zone propagation within the temporal resolution of the method. Conversely, in the heat-treated sample (orange) β decreased throughout, reflecting an increase in the relative proportion of larger microcracks. This change occurred gradually at first and then more sharply once cracks localized along the optimally oriented shear zone, similar to that observed in numerical simulations (Kun et al., 2013) and as inferred from AE magnitude distributions (Sammonds et al., 1992) in dry porous media. The sharp drop in β is a clear precursor to failure, corresponding to propagation of the shear zone through the sample (Figures 5m-5o). This supports our hypothesis that the heat-treated sample exhibits the clear precursors associated with a second-order phase transition.
The evolution of the two-point correlation (fractal) dimension D was very different between the two samples ( Figure 13c). Initially there was a greater degree of clustering in the untreated sample (blue) than the heat-treated one (orange). In both samples D shows a minimum two steps before the onset of localization became apparent in the µCT images, demonstrating the sensitivity of D to localization (see also Bonnet et al., 2001). The degree of clustering at this point, reflected in the value of D, was very similar between the two samples. From this point on, D in the untreated sample increased significantly as microcracks became more distributed (less clustered) along the radial zones. Conversely, in the heat-treated sample D shows increased clustering that was sustained throughout localization. It recovered (implying decreased clustering) to a relatively constant value as the optimally oriented shear zone propagated stably through the sample before accelerating at the point of failure as the shear zone spanned the sample. Thus, D highlights clear differences in the spatial distribution of the microcrack network between the increasingly distributed damage zones in the untreated sample and localization of an asymmetric shear fault in the heat-treated one. The increasingly distributed nature of crack damage in the untreated sample gives no indication of potential failure, while increased clustering due to localization in the heat-treated sample is a clear and early precursor to failure associated with the development of a damage zone optimally oriented to encourage system-sized shear failure. While both samples show precursory changes, only the heat-treated sample has precursors capable of accurately forecasting the point of system-sized catastrophic failure. −p in the heat-treated sample (orange), to a predicted failure stress, σ p , within 3% of the observed failure stress, σ c . This compares with the same model yielding a much poorer prediction in the untreated sample (blue), where σ p = 2.4 σ c . Dotted lines correspond to the exponential model, ξ/l = h exp(qσ), which in the heat-treated sample is less likely and in the untreated sample is more likely than the inverse power law model (the relative likelihood of each model is expressed in the main text and also in Table S5, along with the AICc values). The dash-dot lines show σ c for both samples, and the dashed line shows σ p for the heat-treated sample.

10.1029/2020JB019642
Journal of Geophysical Research: Solid Earth

Discussion
The results presented above reveal key aspects of the evolving nature of compressive failure of brittle rocks through the accumulation of microcracks that spontaneously organize themselves along localized damage zones. Our synchrotron X-ray μCT observations of in situ compressive rock deformation reveal the underlying processes-in particular the nature of the phase transition between intact and failed states in materials with different degrees of starting heterogeneity. Both our post-failure samples contained a localized shear fault, but the preceding accumulation of microcracks was very different between the samples, especially in their spatial distribution and their growth characteristics close to failure. We confirm our hypothesis that, in terms of stress and within the time resolution of our experiments, the transition to failure is abrupt and unpredictable (first-order) in the homogeneous sample, but continuous and predictable (second-order) in the heterogeneous sample.

Microcrack Network Evolution
Prior to failure, our initially crack-free, and therefore more homogeneous, sample accumulated damage in a spalling pattern of localized zones distributed radially around the sample with no preferred strike direction. This damage pattern was completely overprinted during failure, highlighting the drawback of analyzing failed samples retrospectively to gain insights into pre-failure damage accumulation. Pre-failure behavior in this sample resembles strain localization observed from in situ μCT images of deforming mono-minerallic, fine-grained, and uniformly graded (i.e., structurally homogeneous) sand specimens (Desrues et al., 1996). The macroscopic fault localized abruptly at >97% of peak (failure) stress, σ c , as microcracks transitioned from being broadly distributed throughout the sample (albeit along several radially oriented zones) to being organized along an emerging shear zone.
In contrast, our precracked, and therefore more heterogeneous, sample accumulated damage around, and subsequently failed along, a localized shear zone. This behavior resembles the observations of Lockner et al. (1991Lockner et al. ( , 1992 who showed progressive localization of AE along a shear zone in deforming Westerly granite samples from peak stress onward. However, in our experiment the shear zone localized earlier, at 90% of σ c , with a subsequent period of stable crack nucleation and growth along the damage zone during strain hardening prior to dynamic rupture at peak stress. This behavior resembles fault nucleation and propagation from AE in Berea (Lockner et al., 1992) and Clashach (Lennartz-Sassinek et al., 2014;Liakopoulou-Morris et al., 1994) sandstone samples (arguably more heterogeneous than granite samples in terms of their porosity), where a diffuse damage zone appeared and gradually localized around an incipient fault plane prior to σ c .
Our results show that heterogeneity exerts a strong control on the evolution of crack network anisotropy, with homogeneity acting to stabilize the system prior to dynamic failure, generating more complex patterns of strain localization with more isotropic global characteristics, as suggested by Desrues et al. (1996). Under axi-symmetric tri-axial loading conditions, sample homogeneity is a constraint that favors a transversely isotropic spalling pattern until very close to peak stress, whereas the presence of heterogeneity acts to amplify the pre-existing anisotropy with the formation of a shear fault. Radial spalling patterns are rarely observed in studies of AE, potentially due to limits on their location accuracy, where microcracks

10.1029/2020JB019642
Journal of Geophysical Research: Solid Earth occurring along several radially distributed, but localized, damage zones might give the impression of being distributed throughout the sample.
In both of our samples, damage accumulated via the nucleation and subcritical growth of microcracks along localized damage zones. En-echelon and wing-crack arrays formed at different stages in the deformation process in each sample (at initial localization in the untreated sample but only once the optimally oriented shear zone localized in the heat-treated sample) and formed at the same degree of strain (Figures 4f and 5l), implying a degree of strain control. The main direction of individual microcrack growth in the localized zones was along strike rather than down dip (Figures 10a-10c). Models of damage accumulation under tri-axial compression are usually based on AE locations and microstructural observations of post-failure samples, from which it is difficult to quantify the relative proportion of progressive, pre-failure axial to radial microcrack growth. Along strike growth is consistent with our conventional tri-axial compressional stress configuration (σ 1 > σ 2 = σ 3 ), in which it is energetically more favorable for tensile microcracks to open radially against the axes of minimum principal stress and close against the axis of maximum principal stress. Down dip fault propagation occurred instead by the nucleation, growth, and then linkage of an increasing number of small, tensile en-echelon and wing-cracks forming at the fault tip (Figures 6a, 6b, 7, and 9c). This is consistent with previous experimental and modeling work (e.g., Ashby & Hallam, 1986;Ashby & Sammis, 1990;Cho et al., 2007;Horii & Nemat-Nasser, 1985Kranz, 1979;Nemat-Nasser & Horii, 1982;Nemat-Nasser & Obata, 1988;Potyondy & Cundall, 2004;Reches & Lockner, 1994;Rundle & Klein, 1989;Sammis & Ashby, 1986;Tapponnier & Brace, 1976) and recent in situ observations of damage accumulation in strong rocks (Renard et al., 2017(Renard et al., , 2018. We observed significant anisotropy of void strike in the pre-existing porosity in both samples (Figures 6c and  6d), despite visual inspection of thin sections and compressional wave velocities of the same rock type showing only 1% anisotropy in bench-top tests on the original material (Meredith et al., 2005;Meredith, pers comm. 02 December 2019). This indicates that a small velocity anisotropy represents substantial void anisotropy. The pre-existing void anisotropy is more pronounced in the heat-treated sample than in the untreated sample, possibly due to thermal expansion during the heat-treatment acting to close the isolated, mainly round voids in the feldspar microphenocrysts . This may also account for the otherwise counterintuitive smaller overall porosity in the heat-treated sample compared with the untreated one. The application of confining pressure may also have contributed to the porosity difference by acting to close the thermally induced cracks in the heat-treated sample more effectively than the stiffer pores in the untreated sample. In the heterogeneous sample, the preferred strike of the pre-existing porosity corresponds almost exactly to the strike of the emerging fault plane (Figure 6d). There was also significant amplification of the pre-existing anisotropy of the rock fabric (from 33% to 96% just before failure; Table S1). This was not the case in the homogeneous sample, where the degree of anisotropy remained approximately constant throughout the lead-up to failure (Table S1), consistent with the lack of an overall preferred strike in the pre-failure localized zones in this sample.
The results in the previous paragraph prove that the initial microstructure, specifically the orientation and anisotropy of pre-existing porosity dictated the geometry and location of the future (post-failure) fault, particularly in the heat-treated sample. We speculate that this happens via a modification of the local stress field with respect to the principal stress axes. In true tri-axial configurations (σ 1 > σ 2 > σ 3 ), shear wave velocity anisotropy measurements have shown that microcracks in general propagate parallel to σ 2 as they open parallel to σ 3 (Crawford et al., 1995), while polymodal faulting is also often seen (Healy et al., 2015). Thus, although the global stress configuration is axi-symmetric in our case, both heterogeneity and void anisotropy in the microstructure appear to cause the local development of truly tri-axial stresses such that a particular strike is preferred. One possible mechanism for this may be stress rotation around microstructural discontinuities (Faulkner et al., 2006), possibly reflected in our experiments in the rotation of the void ellipsoids with respect to the principal stress axes (Figures 6c and 6d). In this case, the pre-existing network of anisotropic microcracks with a preferred orientation would have generated an emergent, locally dominant true-tri-axial stress field within the body of our heterogeneous sample, even though the confining pressure was isotropic around the vertical (σ 1 ) axis. Conversely, in our homogeneous sample, some complex interplay between local true tri-axial stresses and global axi-symmetry would be required to generate several radially distributed damage zones simultaneously. We speculate that the global axi-symmetry initially counteracts the rotation of internal stresses in this sample, acting to prevent an increase in crack anisotropy and thereby increasing the uniformity of the strike distribution as the experiment progresses. Thus, the relationship between the evolving anisotropy of the microcracks and their preferred orientation is likely to be a controlling factor on the geometry and location of an asymmetric shear fault, on the timing of the formation of this fault and on whether pre-failure damage is localized along this fault or not.
In both our samples, the majority of cracks dip steeply within ±15°of the loading direction, although a few dip less steeply between 15°and 30° (Figures 6a and 6b). This is consistent with the results of post-failure sample analysis in early experimental work (Brace et al., 1966;Hallbauer et al., 1973;Lajtai, 1974). The macroscopic fault in our homogeneous sample dips at a similar angle to the pre-failure microcracks, whereas in our heterogeneous sample it dips less steeply post-failure than it does at peak stress. Although the effective pressure was relatively low (15 MPa), which may promote axial failure over shear, it was consistent across the two experiments. This implies that the differences in fault dip result from an intrinsic microstructural response, whereby the emergent internal friction coefficient decreases during failure in the heterogeneous sample but remains constant in the homogeneous sample, consistent with DEM models (Kun et al., 2018) that show a decreasing coefficient of internal friction with increasing heterogeneity. In both samples, the dip angle increases during quasi-static damage accumulation, increasing earlier in the homogeneous case (during initial localization) than the heterogeneous case (only during localization around the optimally oriented shear zone). In the homogeneous case, the steep dip of the nucleating cracks (Figures 6a and 7a) and the eventual fault plane (10°; Figure 6a) indicates that the internal friction coefficient in this sample is sufficiently high to inhibit microcrack damage by shear mechanisms until immediately before dynamic failure. In the heterogeneous case, the dip, and therefore the internal friction coefficient, increases only during propagation of the shear zone and is particularly pronounced immediately before failure (Figure 6b), while the eventual fault plane dips less steeply (30°; Figure 6b). This indicates that early crack nucleation and failure itself both involve shear mechanisms, but pre-failure shear zone propagation is governed primarily by tensile mechanisms, that is, the accumulation of en-echelon tensile cracks (Figure 7b), with a corresponding increase in the internal friction coefficient. For this reason we have referred to a "damage zone" prior to failure and a "fault plane" afterward.
In our homogeneous sample, increased clustering (Figure 13c; blue circles) occurred at 43% σ c with the onset of localization at 64%. This agrees with observations and models of cracks initiating earlier than the theoretical shear-sliding threshold for more homogeneous low porosity, crystalline rocks (70% σ c ; Hallbauer et al., 1973;Nicksiar & Martin, 2013. The implication is that our more homogeneous sample is weakest in tension and, once a sufficient number of tensile cracks form, a macroscopic shear fracture will naturally develop. We therefore conclude that damage in this sample most likely initiated via the nucleation of pore-emanating (Ashby & Sammis, 1990;Sammis & Ashby, 1986) or force-chain controlled (Cho et al., 2007;Potyondy & Cundall, 2004) tensile microcracks due to the redistribution of stress around equant compressing pores and grains. Conversely, increased clustering in our heterogeneous sample (Figure 13c; orange circles) occurred at 62% σ c with the onset of localization at 72%. This is later than the theoretical shear-sliding threshold for heterogeneous rocks (60% σ c ; Hallbauer et al., 1973;Nicksiar & Martin, 2013. The implication here is that our more heterogeneous sample is weaker in shear than in tension since shear-sliding along preferentially oriented, pre-existing cracks occurred before tensile cracking. We therefore conclude that damage in this sample most likely initiated via the development of tensile "wing-cracks" (Ashby & Hallam, 1986;Ashby & Sammis, 1990;Horii & Nemat-Nasser, 1985Nemat-Nasser & Horii, 1982;Nemat-Nasser & Obata, 1988) at the tips of pre-existing defects due to shear-sliding along those defects. Unfortunately, such shear-sliding would not be visible in our images without significant dilatancy during slip.
In summary, our experimental data confirm that the initial heterogeneity within a rock sample is a key control over how cracks, pores, and grain boundaries interact locally with the applied stress field and imply that the microstructure transitions from being weakest in tension to being weakest in shear as heterogeneity increases.

Scaling, Phase Transition Style and Predictability of Failure Time
Microcrack volume and intercrack length distributions follow power laws throughout the cycle of deformation and failure in both samples, characteristic of the scale-invariant (fractal) nature of natural fault networks (Bonnet et al., 2001;Main et al., 1990) and consistent with the power law microcrack volume distributions observed by Renard et al. (2017Renard et al. ( , 2018. The transition from the TRP to the GR model for the microcrack volume distributions ( Figure S2) in the homogeneous sample emulates changes in the organization of earthquake size distributions following the occurrence of extreme or very large earthquakes (Bell et al., 2013a). Close to failure the void volume distribution shows a bump at large volumes, indicative of a supercritical state with an elevated probability of occurrence of large events (Main, 1996), sometimes known as "dragon kings" (Sornette, 2009). We have demonstrated that the parameters of these distributions are more sensitive to heterogeneity than porosity alone, consistent with the findings of Vasseur et al. (2017) and Kun et al. (2018). In combination with μCT observations of fault formation, the evolution of these parameters provides a microstructural explanation for the variation in the systematic prediction error for the failure time based on acoustic emissions (Vasseur et al., 2015).
However, the systematic change in the mean void aspect ratios during crack growth may indicate that the scaling of crack growth is self-affine (i.e., exhibits scale-invariance in length with different exponents for individual growth axes, leading to a variable aspect ratio) rather than self-similar (the same scaling exponent for all growth axes, with a constant aspect ratio). This is consistent with observations of fracture surface geometries in rocks (Schmittbuhl et al., 1995) and other materials (Bouchaud et al., 1990;Mandelbrot et al., 1984;Russ, 1994;Schmittbuhl & Maloy, 1997; see also Bouchaud, 1997, for a review), which are well-described by self-affine fractals. These studies have shown that scaling along the aperture axis is systematically smaller than along the mean crack plane, with the systematic (Hurst) exponent defining the fracture roughness (Bouchaud, 1997;Weiss, 2001). Our observation that almost no growth at all occurs along the aperture axis supports the conjecture that the aperture direction is not physically equivalent to the mean crack plane (Schmittbuhl et al., 1995). Our results indicate that scaling along the strike and dip axes may also systematically differ from each other. This contradicts the notion of strict self-similarity in the mean crack plane (Schmittbuhl et al., 1995) and implies that the strike and dip directions are not physically equivalent either. Further work is required to quantify the scaling anisotropy for crack growth in our experiments and to test these hypotheses. Since crack surfaces in crystalline materials require heterogeneities, such as grain boundaries and dislocations that pin the propagating crack front, in order to develop self-affine roughness (Bouchaud, 1997;Schmittbuhl & Maloy, 1997;Weiss, 2001), we expect that scaling exponents for the heterogeneous sample may be more anisotropic than for the homogeneous sample.
In the heterogeneous (heat-treated) case, we find evidence for a continuous (second-order) phase transition in the inverse power law acceleration to failure of ξ with respect to stress (Figure 12; solid orange line), with failure occurring near the asymptote, together with clear precursors in β and D. The rapid decrease in β corresponds to the formation of a localized damage zone optimally oriented for macroscopic shear failure, occurring when the microcrack network self-organizes. This provides a clear precursor to sample failure related to a distinct physical process, that is, the emergent inverse power law acceleration in ξ. The asymptote defines a predictable failure time defined by a smooth transition to an infinite ξ at the sample-scale (Figure 1; orange line). The early and sustained decrease in D in 3-D is a key precursory indicator of localization, while its recovery is associated with shear zone propagation in 2-D, as anticipated by the model of Main (1992). This provides another clear precursor to failure. Such behavior agrees with statistical physics models of rupture as a critical, second-order phenomenon (Girard et al., 2010;Kun et al., 2013). Thus, taken together, these variables show that damage localization along a zone optimally orientated for macroscopic shear failure is the physical process that defines whether the phase transition from an intact to a failed state is second-order, and therefore predictable, with reliable precursors to failure.
In the homogeneous (untreated) case, we find evidence for an abrupt or discontinuous (first-order) phase transition, with an unsuccessful forecast of the failure stress, and a preference for an exponential model for the evolution of the correlation length, ξ, with respect to stress. Furthermore, there is very little evidence for reliable precursors in either the microcrack volume exponent, β, or the two-point correlation dimension, D, and the bump in the void size distribution at large volumes is reminiscent of a first-order phase transition (Ceva & Perazzo, 1993;Lominitz-Adler et al., 1992). Approaching failure we see small fluctuations in β, ξ, and D that may indicate impending failure as they are associated with formation of the additional damage zones and subsequent microstructural instability due to crack nucleation close to failure. However, using these parameters as precursors may lead to false alarms since they are not associated with the eventual fault plane. The exponential increase in ξ (implying that local correlations dominate) is unusual and generally associated with the critical regime during phase transitions across surfaces (Kosterlitz, 1974;Kosterlitz & Thouless, 1973), such as during large-scale faceting at the surfaces of growing crystals (Nozieres, 1992). Its

10.1029/2020JB019642
Journal of Geophysical Research: Solid Earth stabilization to a finite value shortly followed by abrupt failure is characteristic of a first-order phase transition (Figure 1; green line). In numerical models of fault growth, an exponential distribution of fault lengths is associated with crack nucleation, whereas a power law distribution emerges with nucleation plus crack growth and coalescence (Cowie et al., 1995). Hence, the origin of this response in our rock volume may be explained by our observation that crack nucleation is the dominant damage process in the homogeneous sample while crack growth becomes increasingly important closer to failure in the heterogeneous sample (section 3.3.1). This behavior corresponds to the existence of a metastable state of crack nucleation at a system-sized ξ during a first-order transition, when the system is vulnerable to the influence of sufficiently large perturbations (subcritical bifurcation) (Sornette, 2006). This vulnerability and the resulting discontinuity may be the reason for an unpredictable failure time (Vasseur et al., 2015).
An estimate for the correlation length exponent (1.15) for Carrara marble (Kandula et al., 2019) falls almost exactly halfway between the exponents for Ailsa Craig microgranite found here (0.65 for the heterogeneous sample and 1.75 for the homogeneous sample). However, the Carrara marble exponent was estimated by assuming the failure stress a priori, so it is not directly comparable with our results. It is therefore not possible to confirm whether an inverse power law would successfully forecast the failure stress in real time and/or whether a different model would be more likely. Nevertheless, the nature of Carrara marble may place it halfway between our two end members. It is chemically pure, composed of 99% annealed calcite crystals (Alber & Hauptfleisch, 1999), with a homogeneous microstructure (Oesterling, 2004), a very low permeability (10 −19 m 2 ) and only 0.2-0.5% connected porosity (Alber & Hauptfleisch, 1999;Bandini et al., 2012;Cartwright-Taylor, 2015;Zhang et al., 1994). However, studies have shown the presence of microdiscontinuities within grains, including twin lamellae (Bandini et al., 2012;Cartwright-Taylor, 2015;Ramez & Murrell, 1964) and a high density of dislocations (Fredrich et al., 1989), while its isotropic texture consists of both well-locked (xenoblastic) and more mobile (granoblastic) grain boundaries (Bandini et al., 2012;Cartwright-Taylor, 2015). These factors indicate a complex history of both static and dynamic recrystallization (Molli & Heilbronner, 1999;Oesterling, 2004) and introduce a degree of heterogeneity that may be intermediate between our two samples.
In both samples, the critical value of ξ is 200 μm, marking the longest crack supported by the sample volume without a runaway instability developing. Significantly, this falls just short of the mean grain size of the groundmass (250 μm). That is, catastrophic failure occurs when whole grains break. This confirms the findings of Vasseur et al. (2017) from acoustic emissions (AE) data that the grain size (inter-particle distance) is a better metric for the characteristic void dimension at failure than the distance between pores (inter-void distance).
Our observations highlight the strong dependence of the degree of predictability on material properties that may be unknown in a field application, as well as the importance of analyzing several independent parameters for identifying the type of phase transition and predicting the point of failure (Lei & Satoh, 2007). They may also explain why, when looking at long time-series of field-scale seismicity or deformation, clear and reliable precursors to failure are detected only in some cases and preferentially in application to forecasting of landslides and volcanic eruptions. In other cases, notably in forecasting of individual large earthquakes, fluctuations related to instability may be present but may not be statistically significant enough to be detectable as precursors. In both samples, D shows increased clustering earlier than localization is visually apparent in the μCT images and therefore may provide useful information about the impending onset of damage localization for a variety of applications and settings. Finally, the relatively high strain rates analyzed here may not be representative of the evolution of precursors at lower strain rates. For example Ojala et al. (2004) showed that the acceleration to failure in AE rate asymptotically approaches the behavior expected of a single Griffith crack (Figure 1) as strain rate is decreased in laboratory compression tests on porous sandstones. Nevertheless, we have confirmed that heterogeneity plays a significant role in determining the style of evolution of the population of microcracks and hence the predictability of the system-scale failure time.

Suggestions for Future Work
This discussion has highlighted some outstanding research questions to be addressed in future work. The most notable of these are as follows: (i) Why do previously obtained degrees of anisotropy inferred from acoustic measurement differ markedly from our newly obtained structural ones? (ii) How does crack growth scale (in terms of the ellipsoid radii), and is it self-affine? (iii) Does the predominant local failure mechanism change from tensile to shear as system-sized failure approaches, as seen in the AE data of Graham et al. (2010)? (iv) Does this transition occur earlier in more heterogeneous materials? Given tensile fractures are easier to see in imaging void space, the latter two questions would benefit from digital volume correlation techniques that are the subject of ongoing work and can detect local changes in shear and volumetric strain.

Conclusions
Our in situ time-resolved X-ray μCT images of very small samples of deforming granite show that the heterogeneity of the starting material exerts a strong control on the evolution of the statistical properties of crack size and spatial distribution during fracture network evolution. The accumulating microcracks have power law frequency-volume and intercrack length distributions over a finite scaling range, irrespective of the degree of starting heterogeneity, with well-determined scaling exponents β (the size exponent) and D (the correlation dimension). The inferred correlation length ξ increases exponentially with respect to stress in the homogeneous case, with sudden-onset, unpredictable failure, analogous to the behavior expected for a first-order (discontinuous) phase transition. In contrast, the heterogeneous sample shows an inverse power law acceleration to a predictable failure point at the asymptote, diagnostic of a second-order (continuous) phase transition (Equation 1 and Figure 1). The second-order transition is linked to the distinct physical process of quasi-static, asymmetric accumulation of damage within an optimally oriented zone increasingly localized around the eventual fault plane, with associated reliable precursors to failure in the evolution of β and D. The correlation dimension is a key early indicator of localization on such a shear zone for the heterogeneous sample. This is not observed within the time resolution of our observations for the homogeneous sample, where the precursory damage takes the form of more radially isotropic zones of spalling, and shear localization must occur very close to the point of dynamic failure itself to explain the post-failure observation of a shear fault.
Crack nucleation dominates the frequency-size statistics in the homogeneous case and crack growth in the heterogeneous case. In both cases, the transition to localized damage occurs by a combination of nucleation and growth. The timing of the onset of crack coalescence defines the order of the phase transition and hence the predictability of the failure time. Nevertheless, catastrophic failure occurs in both cases as the correlation length approaches the grain size, which in turn controls the failure of local bridges between aligned en-echelon and wing-cracks in the shear damage zone in the heterogeneous sample. The initial rock microstructure, specifically the anisotropy of pre-existing porosity, dictates the geometry and orientation of the emergent fault plane, independent of starting heterogeneity. This reflects the strong control of starting microstructure on the rock's internal stress state, despite the axi-symmetric external loading conditions and the very low anisotropy (1%) inferred from acoustic velocity measurement.

Data Availability Statement
The data supporting our conclusions can be found in main text, in the Supporting Information, and in the μCT data sets held at the NERC repository (http://data.ceda.ac.uk/ngdc/R001693-1). When using the μCT data sets, please cite Cartwright-Taylor et al. (2020) and see https://doi.org/10.5285/0dc00069-8da8-474a-8993-b63ef5c25fb8 for the metadata.