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):
- 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.
- 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.
- 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).