Discussion

The results presented above reveal key aspects of the evolving nature of compressive failure of brittle rocks through the accumulation of micro-cracks that spontaneously organize themselves along localized damage zones. Our synchrotron x-ray micro-tomography (μ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 micro-cracks 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, \(\sigma_{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 pre-cracked, 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. (1991; 1992) who showed progressive localization of AE along a shear zone in deforming Westerly granite samples from peak stress onwards. However, in our experiment the shear zone localized earlier, at 90% of\(\sigma_{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 (Liakopoulou-Morris, et al. 1994; Lennart-Sassinek et al., 2014) 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\(\sigma_{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 triaxial 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 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 sub-critical growth of micro-cracks 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 (Figure 4F and Figure 5L), implying a degree of strain control. The main direction of individual micro-crack growth in the localized zones was along strike rather than down dip (Figure 10a,b,c). 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 micro-crack growth. Along-strike growth is consistent with our conventional tri-axial compressional stress configuration (\(\sigma_{1}>\ \sigma_{2}=\sigma_{3}\)), in which it is energetically more favorable for tensile micro-cracks 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 (Figure 6a,b, Figure 7 and Figure 9c). This is consistent with previous experimental and modelling work (e.g., Tapponnier and Brace, 1976; Kranz, 1979; Nemat-Nasser and Horii, 1982; Horii and Nemat-Nasser, 1985; 1986; Sammis and Ashby, 1986; Ashby and Hallam, 1986; Nemat-Nasser and Obata, 1988; Rundle and Klein, 1989; Ashby and Sammis, 1990; Reches and Lockner, 1994; Potyondy and Cundall, 2004; Cho et al., 2007) and recentin-situ observations of damage accumulation in strong rocks (Renard et al., 2017; 2018).
We observed significant anisotropy of void strike in the pre-existing porosity in both samples (Figure 6c,d), 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.). 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 micro-phenocrysts (Meredith et al., 2012). This may also account for the otherwise counter-intuitive 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 in SI). 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 in SI), 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 (\(\sigma_{1}>\ \sigma_{2}>\sigma_{3}\)), shear wave velocity anisotropy measurements have shown that micro-cracks in general propagate parallel to \(\sigma_{2}\) as they open parallel to\(\sigma_{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 (Figure 6c,d). In this case, the pre-existing network of anisotropic micro-cracks with a preferred orientation would have generated an emergent, locally dominant true-triaxial stress field within the body of our heterogeneous sample, even though the confining pressure was isotropic around the vertical (\(\sigma_{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 micro-cracks 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° (Figure 6a,b). 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 micro-cracks, 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 (Figure 6a and Figure 7a) and the eventual fault plane (10°; Figure 6a) indicates that the internal friction coefficient in this sample is sufficiently high to inhibit micro-crack 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 shear zone propagation is governed primarily by tensile mechanisms, i.e., 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’ afterwards.
In our homogeneous sample, increased clustering (Figure 13c; blue circles) occurred at 43% \(\sigma_{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% \(\sigma_{c}\); Hallbauer et al., 1973; Nicksiar and Martin, 2013; 2014). 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 (Sammis and Ashby, 1986; Ashby and Sammis, 1990) or force-chain controlled (Potyondy and Cundall, 2004; Cho et al., 2007) tensile micro-cracks due to the re-distribution of stress around equant compressing pores and grains. Conversely, increased clustering in our heterogeneous sample (Figure 13c; orange circles) occurred at 62%\(\sigma_{c}\) with the onset of localization at 72%. This is later than the theoretical shear-sliding threshold for heterogeneous rocks (60% \(\sigma_{c}\); Hallbauer et al., 1973; Nicksiar and Martin, 2013; 2014). 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’ (Nemat-Nasser and Hori, 1982; Horii and Nemat-Nasser, 1985; 1986; Ashby and Hallam, 1986; Nemat-Nasser and Obata, 1988; Ashby and Sammis, 1990) 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

Micro-crack volume and inter-crack 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 (Main et al., 1990; Bonnet et al., 2001) and consistent with the power-law microcrack volume distributions observed by Renard et al. (2017; 2018). The transition from the TRP to the GR model for the micro-crack volume distributions (Figure S2 in SI) 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 (Mandlebrot et al., 1984; Bouchaud et al., 1990; Russ, 1994; Schmittbuhl and 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 (Schmittbuhl and Maloy, 1997; Bouchaud, 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 \(\xi\) with respect to stress (Figure 12; solid orange line), with failure occurring near the asymptote, together with clear precursors in \(\beta\) and \(D\). The rapid decrease in\(\beta\) 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, i.e. the emergent inverse power-law acceleration in \(\xi\). The asymptote defines a predictable failure time defined by a smooth transition to an infinite\(\xi\) at the sample-scale (Figure 1; orange line). The early and sustained decrease in \(D\) in 3D is a key precursory indicator of localization, while its recovery is associated with shear zone propagation in 2D, 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, \(\xi\), with respect to stress. Furthermore, there is very little evidence for reliable precursors in either the micro-crack volume exponent, \(\beta\), or the two-point fractal dimension, \(D\), and the bump in the void size distribution at large volumes is reminiscent of a first-order phase transition (Lominitz-Adler et al., 1992; Ceva and Perazzo, 1993). Approaching failure we see small fluctuations in \(\beta\), \(\xi\) 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 \(\xi\) (implying that local correlations dominate) is unusual and generally associated with the critical regime during phase transitions across surfaces (Kosterlitz and Thouless, 1973; Kosterlitz, 1974), such as during large-scale faceting at the surfaces of growing crystals (Nozières, 1992). Its 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.2.1). This behavior corresponds to the existence of a metastable state of crack nucleation at a system-sized \(\xi\) 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 stressa 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 and Hauptfleisch, 1999), with a homogeneous microstructure (Oesterling, 2004), a very low permeability (10-19m2) and only 0.2-0.5% connected porosity (Zhang, 1994; Alber and Hauptfleisch, 1999; Bandini et al., 2012; Cartwright-Taylor et al., 2015). However, studies have shown the presence of micro-discontinuities within grains, including twin lamellae (Ramez and Murrell, 1964; Bandini et al., 2012; Cartwright-Taylor et al., 2015) 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 et al., 2015). These factors indicate a complex history of both static and dynamic recrystallization (Molli and 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 \(\xi\) 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 and 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 behaviour 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 micro-cracks, 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.