Figure 3: (a ) Segmented fracture network (blue) within a greyscale 3D μCT sub-volume at peak stress in the heat-treated sample. (b ) Example of a 2D grey-scale slice, also with segmented fractures highlighted in blue. Most of the very narrow fractures remain undetected due to the similarity of the grey-scale values with the surrounding material.

Analysis of the segmented microcrack network

Porosity and the number of fractures

Initially of interest is the evolution of porosity, \(\varphi\)(including both the pre-existing pores and the developing micro-cracks) and the total number of voids, \(N\), in each 3D sub-volume 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 AvizoTM. Porosity was obtained from the total void volume divided by the analyzed sub-volume: \(\varphi=\ \frac{(\sum{V_{i})}}{(\pi R^{2}l)}\), where \(V_{i}\) is volume of each crack, \(R\) is the sample radius and\(l\) is the length of the analyzed sub-volume. To determine the most likely empirical relationship with strain for both \(\varphi\) and\(N\), we found the parameters for several possible models (quadratic, exponential and simple power-law) using non-linear least-squares regression and then used the corrected Akaike Information Criterion,AICc (Hurvich and Tsai, 1989) to test these competing models objectively (see Text S2 in the SI for a full description of the calculations).

3D micro-crack 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 micro-crack). Each ellipsoid was calculated independently from the crack’s 3D moments of inertia (Text S1b in SI), also using the Label Analysis module in AvizoTM. 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=\sqrt{5\left|\text{eigenvalue}\right|}\) 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, \(\xi\) (linear dimension of the largest void), the void size exponent, \(\beta\), 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; Sykes and Jaumé, 1990; Bufe and Varnes, 1993; Sornette and Sammis, 1995; Turcotte, 1997; Ouillion and Sornette, 2000; Zöller et al., 2001; Kagan, 2002; Sammis and Sornette, 2002; Tyupkin and Giovambattista, 2005). In rock deformation studies (e.g., Moura et al., 2005; Lei and Satoh, 2007), these parameters have been similarly inferred from the amplitudes and locations of acoustic emissions (AE). In particular, the inter-crack 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 \(\xi\), \(\beta\) and \(D\) directly from the evolving population of micro-cracks 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 in SI) directly from the Label Analysis module in AvizoTM, and then computed Euclidean lengths, \(L_{i}\), between centroids.
We obtained maximum likelihood estimates of the void size exponent,\(\beta\), from the frequency-volume distribution in each µCT sub-volume. 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 in SI):
  1. GR: the Pareto distribution (a pure power law, equivalent to the Gutenberg-Richter distribution) with cumulative complementary (survivor) function \(\Phi\left(V\right)={(V_{t}/V_{i})}^{\beta}\)for \(V_{t}<V_{i}<\ V_{\max}\), where \(V_{i}\) is volume of each individual void and \(V_{\max}\) is the upper bound (maximum) void volume in the distribution.
  2. TRP: the truncated Pareto distribution (similar to the GR but showing a power-law taper in the tail towards \(V_{\max}\)), with\(\Phi\left(V\right)={[(V_{t}/V_{i})}^{\beta}-\ \left(\frac{V_{t}}{V_{c}}\right)^{\beta}]/[1-\ \left(\frac{V_{t}}{V_{c}}\right)^{\beta}]\)for \(V_{t}<V_{i}<\ V_{c}\), where \(V_{c}\) is the tapering corner volume of the distribution.
  3. TAP: the tapered Pareto distribution (equivalent to the modified Gutenberg-Richter relation which shows an exponential taper in the tail towards \(V_{\max}\)), with\(\Phi\left(V\right)={(V_{t}/V_{i})}^{\beta}\exp{[(V_{t}-V_{i})/V_{c}]}\)for \(V_{t}<V_{i}<\ \infty\).
We defined a correlation length, \(\xi=\sqrt[3]{V_{\max}}\) for the pure power-law model or \(\sqrt[3]{V_{c}}\) 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 (Rydelek and Sacks, 1989; Woessner and Wiemer, 2005; Mignan and Woessner, 2012), and is equivalent to the threshold of completeness in seismicity data. We obtained \(V_{t}\) from the maximum curvature method (Roberts et al., 2015), i.e., 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 SI). In both samples, the number of voids with \(V\geq V_{t}\) exceeded 200, which is the minimum catalogue size required for reliable estimation of \(\beta\) (Roberts et al., 2015). We used a modified Bayesian Information Criterion (BIC) (Leonard and Hsu, 1999; Bell et al., 2013a) to find the most appropriate model (see Text S3 for full details of the calculation, and Figure S2 in SI for the full results) thus obtained the most likely values of \(\beta\) and \(\xi\) for each distribution. TheBIC 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, i.e., whether or not\(\xi\) followed an inverse power-law acceleration to failure, we fit the data by non-linear 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 non-linear regression. The point of failure,\(x_{p}\), is defined by a mathematical singularity as\(y\rightarrow\infty\). 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., Voight, 1988; Cornelius and Voight, 1994; Utsu et al., 1995; Kilburn and Voight, 1998; Kilburn 2003; Smith et al., 2009; Bell et al., 2013b; Vasseur et al., 2015; 2017). 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, \(\sigma_{p}\), is what would be available in real time, rather than the observed failure stress, \(\sigma_{c}\). By estimating \(\sigma_{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 in SI for details). We also tested an exponential model\(y=h\ \exp(qx)\) as an alternate hypothesis. We cannot use a simple criterion such as r2 alone to determine the relative goodness of fit because the competing hypotheses have different degrees of freedom, so we used the corrected Akaike Information Criterion (AICc – see Text S2 in the Supporting Information for details of the calculations). It is based on the residual sum of squares, and is considered more robust than r2 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 fractal dimension, \(D\), from the relation\({P(L}_{i})\ \sim\ {L_{i}}^{D-1}\) ( Turcotte, 1997), where\({P(L}_{i})\) is the incremental probability distribution of the inter-void lengths, \(L_{i}\) (see Text S4 in the SI for more detail). The exponent, \(D_{\text{inc}}=D-1\), of \({P(L}_{i})\) in the identified power-law region, \(30<L_{i}<1350\) μm, was obtained from a linear regression in log-log space (Figure S3 in the SI). \(D\) is then \(D_{\text{inc}}+1\). If \(D\ <\ 3\), voids are clustered spatially but as \(D\ \rightarrow\ 3\) they become volume-filling, and therefore less clustered (Robertson et al., 1995).