\documentclass[10pt]{article}

\usepackage{fullpage}
\usepackage{setspace}
\usepackage{parskip}
\usepackage{titlesec}
\usepackage[section]{placeins}
\usepackage{xcolor}
\usepackage{breakcites}
\usepackage{lineno}
\usepackage{hyphenat}





\PassOptionsToPackage{hyphens}{url}
\usepackage[colorlinks = true,
            linkcolor = blue,
            urlcolor  = blue,
            citecolor = blue,
            anchorcolor = blue]{hyperref}
\usepackage{etoolbox}
\makeatletter
\patchcmd\@combinedblfloats{\box\@outputbox}{\unvbox\@outputbox}{}{%
  \errmessage{\noexpand\@combinedblfloats could not be patched}%
}%
\makeatother


\usepackage{natbib}




\renewenvironment{abstract}
  {{\bfseries\noindent{\abstractname}\par\nobreak}\footnotesize}
  {\bigskip}

\titlespacing{\section}{0pt}{*3}{*1}
\titlespacing{\subsection}{0pt}{*2}{*0.5}
\titlespacing{\subsubsection}{0pt}{*1.5}{0pt}


\usepackage{authblk}


\usepackage{graphicx}
\usepackage[space]{grffile}
\usepackage{latexsym}
\usepackage{textcomp}
\usepackage{longtable}
\usepackage{tabulary}
\usepackage{booktabs,array,multirow}
\usepackage{amsfonts,amsmath,amssymb}
\providecommand\citet{\cite}
\providecommand\citep{\cite}
\providecommand\citealt{\cite}
% You can conditionalize code for latexml or normal latex using this.
\newif\iflatexml\latexmlfalse
\providecommand{\tightlist}{\setlength{\itemsep}{0pt}\setlength{\parskip}{0pt}}%

\AtBeginDocument{\DeclareGraphicsExtensions{.pdf,.PDF,.eps,.EPS,.png,.PNG,.tif,.TIF,.jpg,.JPG,.jpeg,.JPEG}}

\usepackage[utf8]{inputenc}
\usepackage[ngerman,greek,english]{babel}



\usepackage{float}








\begin{document}

\title{Catastrophic failure: how and when? Insights from 4D in-situ x-ray
micro-tomography}



\author[1,1,1]{Alexis Cartwright-Taylor}%
\author[1,1,1]{Ian G Main}%
\author[1,1,1]{Ian B Butler}%
\author[1,1,1]{Florian Fusseis}%
\author[1,1,1]{Michael Flynn}%
\author[2,2,2]{Andrew King}%
\affil[1]{University of Edinburgh}%
\affil[2]{SOLEIL Synchrotron, Gif-sur-Yvette, France}%


\vspace{-1em}



  
  \date{May 25, 2024}


\begingroup
\let\center\flushleft
\let\endcenter\endflushleft
\maketitle
\endgroup





\selectlanguage{english}
\begin{abstract}
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 micro-cracks that spontaneously
self-organize along localized damage zones under compressive stress.
Here we present x-ray micro-tomography 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 micro-granite; known for being virtually crack-free) on crack
network evolution and localization. To control for heterogeneity, we
introduced a random nano-scale 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.%
\end{abstract}%



\sloppy


Alexis Cartwright-Taylor\textsuperscript{1}, Ian G.
Main\textsuperscript{1}, Ian B. Butler\textsuperscript{1}, Florian
Fusseis\textsuperscript{1}, Michael Flynn\textsuperscript{1} and Andrew
King\textsuperscript{2}

\textsuperscript{1}School of Geosciences, University of Edinburgh, James
Hutton Road, Edinburgh EH9 3FE, UK. \textsuperscript{2}SOLEIL
Synchrotron, L'Orme des Merisiers Saint Aubin, Gif-sur-Yvette, France.

Corresponding author: Alexis Cartwright-Taylor
alexis.cartwright-taylor@ed.ac.uk

\section*{Key Points:}

{\label{key-points}}

\begin{itemize}
\tightlist
\item
  Micron-scale crack network evolution and localization is
  imaged\emph{in-situ} in deforming granite samples with controlled
  initial heterogeneity
\item
  Heterogeneity affects order and predictability of the phase transition
  to failure, which occurs as correlation length approaches grain size
\item
  Correlation dimension of the crack population is a key indicator of
  localization, and anisotropy in starting porosity dictates fault
  strike
\end{itemize}

\section*{Abstract}

{\label{abstract}}

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 micro-cracks that spontaneously
self-organize along localized damage zones under compressive stress.
Here we present x-ray micro-tomography observations that elucidate
the\emph{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 micro-granite; known for being virtually crack-free) on crack
network evolution and localization. To control for heterogeneity, we
introduced a random nano-scale 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.

\section*{Plain Language Summary}

{\label{plain-language-summary}}

When rocks deform, tiny cracks appear, increasing in size and number
until the rock breaks completely, often along a narrow plane of weakness
where cracks have spontaneously aligned. Sometimes, when the
microstructure is complicated, cracking accelerates quickly in a
predictable way, giving a good indication of when the rock will break.
In other cases, when the microstructure is more uniform, cracking
accelerates more slowly and the rock breaks suddenly and early. To
understand why failure is predictable in some cases but not others -- a
major problem in managing risk from material failure (e.g., earthquakes)
-- we used x-ray imaging to see how cracks form and interact with each
other inside deforming rocks. We found that predictable behavior only
arose when cracks aligned themselves asymmetrically. The orientation of
this damage zone was governed by the rock's pre-existing microstructure.
We also found distinct changes in crack size and spatial arrangement
during alignment, indicating that the rock was approaching failure.
However, when cracks did not align asymmetrically, similar changes were
not observed and failure was not predictable. Our results are important
because they help explain why reliable indicators of catastrophic
failure are not always observed and why forecasting may only work in
certain cases.

\section*{List of Symbols (in order of
appearance)}

{\label{list-of-symbols-in-order-of-appearance}}\selectlanguage{english}
\begin{longtable}[]{@{}ll@{}}
\toprule
\textbf{variable} & \textbf{symbol}\tabularnewline
\midrule
\endhead
correlation length (typical linear dimension of largest void) &
\(\xi\)\tabularnewline
time & \(t\)\tabularnewline
critical (failure) time & \(t_{c}\)\tabularnewline
axial sample strain & \(\epsilon\)\tabularnewline
critical (failure) strain & \(\epsilon_{c}\)\tabularnewline
strain rate & \(\dot{\epsilon}\)\tabularnewline
differential stress & \(\sigma\)\tabularnewline
axial stress & \(\sigma_{1}\)\tabularnewline
confining pressure (radially-symmetric) &
\(\sigma_{3}\)\tabularnewline
sample length & \(L\)\tabularnewline
critical (failure) stress & \(\sigma_{c}\)\tabularnewline
sample co-ordinate system & \(x,y,z\)\tabularnewline
dip & \(\phi\)\tabularnewline
strike & \(\theta\)\tabularnewline
porosity & \(\varphi\)\tabularnewline
number of segmented voids & \(N\)\tabularnewline
individual void & \(i\)\tabularnewline
volume of an individual void (pore or micro-crack) &
\(V_{i}\)\tabularnewline
volume of the largest void & \(V_{\max}\)\tabularnewline
sample radius & \(R\)\tabularnewline
length of analyzed sub-volume & \(l\)\tabularnewline
volume of analyzed sub-volume & \(V_{\text{subvol}}\)\tabularnewline
corrected Akaike Information Criterion &
AICc\tabularnewline
first-order moment of inertia (void centroid) &
\(M\)\tabularnewline
second-order moment of inertia (covariance matrix) &
\(I\)\tabularnewline
void ellipsoid radius & \(r\)\tabularnewline
mean void ellipsoid radius & \(\overset{\overline{}}{r}\)\tabularnewline
power-law exponent of void volume distribution &
\(\beta\)\tabularnewline
Euclidean distance between two voids (inter-void distance) &
\(L_{i}\)\tabularnewline
two-point correlation (fractal) dimension (power-law exponent of
inter-void distance distribution) & \(D\)\tabularnewline
cumulative complementary (survivor) function of frequency-volume
distribution & \(\Phi\left(V\right)\)\tabularnewline
completeness threshold of frequency-volume distribution &
\(V_{t}\)\tabularnewline
tapering corner threshold of frequency-volume distribution &
\(V_{c}\)\tabularnewline
modified Bayesian Information Criterion &
BIC\tabularnewline
incremental probability distribution of inter-void lengths &
\({P(L}_{i})\)\tabularnewline
exponent of incremental inter-void length distribution &
\(D_{\text{inc}}\)\tabularnewline
\bottomrule
\end{longtable}

\section*{Introduction}

{\label{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 and 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 micro-fracturing 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 (Main et al.,
1990; Main, 1996; Bonnet et al., 2001), and in the way their size
relates to the energy released during rupture (Abercrombie, 1995;
Goodfellow and 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 and 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 sub-critically, i.e., without causing
system-sized failure (Lockner et al., 1991; Lei et al., 2000), until the
accumulating cracks self-organize along an asymmetric, localized damage
zone. System-scale failure then occurs when nucleating micro-cracks 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 (Sammonds et al., 1992; Moura et al., 2005; 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 and 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 (Gu\selectlanguage{ngerman}éguen and Schubnel, 2003; Alava et al.,
2008). 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, e.g., temperature (Stanley, 1971,
Fig. 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 towards the
critical point (Bruce and 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
and Sornette, 2002). If this occurs in the lead up to macroscopic
failure, then the failure time can be forecast accurately (Vasseur et
al., 2015; 2017). Inverse power-law acceleration to a well-defined
failure time has also been seen in the evolving microstructure
(micro-crack porosity and the volume of the largest micro-fracture) of
crystalline rocks undergoing brittle deformation (Renard et al., 2017;
2018).

However, the evolution of damage does not always allow a fit to a model
with a well-defined failure time. In structurally homogeneous materials,
there is no emergent, smooth power-law acceleration to failure, as shown
experimentally by Vasseur et al. (2015; 2017) for a range of rock types
and material analogues. In the extreme case of a single flaw in an
otherwise uniform starting material, there is no precursor, and
catastrophic system-sized failure occurs suddenly when the flaw
propagates dynamically at a maximum in the 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; 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, i.e.,
cracking which does not fulfil the conditions for sustained propagation,
shows only an exponential increase in the event rate time-to-failure
behavior (Vasseur et al., 2015; 2017). Failure occurs suddenly and
early; often much earlier than expected from an exponential asymptote
(Vasseur et al., 2015; 2017). This behavior is also characteristic of an
abrupt first-order transition, with the exponential behavior reflecting
local correlations (Stanley, 1971; Sethna, 2006).

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, \(\xi\), 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 \(\xi\) in the
critical state. In this case,\(\xi\) 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
and 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
\(\xi\) 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
\(\xi\) 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., 2015; 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 \(\xi\) 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 and Wallace, 1989):\selectlanguage{english}
\begin{longtable}[]{@{}ll@{}}
\toprule
\(\xi\sim{(T_{c}-T)}^{-\nu}\) & (1)\tabularnewline
\bottomrule
\end{longtable}

This is the type of behavior associated with progressively more
heterogeneous materials undergoing brittle deformation (Vasseur et al.,
2015; 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 micro-cracks 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 micro-crack 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 micro-granite: one being an as-received control
(nominally crack-free), and the other containing a pre-existing
nano-scale 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; 2012). They
are an extreme end member of lowest crack density in natural rocks.
Through the analysis of 4D, \emph{in-situ} synchrotron x-ray
micro-tomography (\selectlanguage{greek}μ\selectlanguage{english}CT) images of the two samples undergoing tri-axial
deformation (see Cartwright-Taylor et al. (2020) for access to the
dataset), 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 pre-cracked sample. \emph{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 micro-fractures.
Using a combination of visual inspection of the \selectlanguage{greek}μ\selectlanguage{english}CT images, geometrical
analysis of the evolving crack network, and techniques used in
statistical seismology, we show that the micro-crack network evolution
varies significantly between the two samples. The degree of starting
heterogeneity controls (i) the evolving spatial clustering and
anisotropy of the micro-cracks, 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 pre-cracked 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,
the power-law scaling of the micro-crack volume and inter-crack length
distributions, and some crack growth characteristics, 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.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image1/image1}

\end{center}
\end{figure}

\textbf{Figure 1:} Schematic graph showing the expected behavior of the
correlation length \(\xi\), normalized by its initial value
\(\xi_{0}\), as strain \(\epsilon\)increases towards
failure at critical strain \(\epsilon_{c}\) for the three phase
transition cases described in the text. The Griffith crack (blue line)
and the exponential (orange line) are both examples of first-order
transitions, exhibiting abrupt failure evident as the discontinuity in
\(\xi\) at \(\epsilon_{c}\). The inverse power-law (green
line) shows the smooth, continuous evolution to failure characteristic
of a second-order transition.

\section*{Materials and Methods}

{\label{materials-and-methods}}

\subsection*{Ailsa Craig micro-granite and thermal
stressing}

{\label{ailsa-craig-micro-granite-and-thermal-stressing}}

We used Ailsa Craig micro-granite (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
(Meredith et al., 2012). As received from the quarry, no pre-existing
microcracks are detectable either by optical or scanning electron
microscopy (Meredith et al., 2005; 2012). Porphyritic in texture with a
groundmass of mean grain size 250 \selectlanguage{greek}μ\selectlanguage{english}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
(\textless{}\textless{} 1 \%) and permeability (1.5 x
10\textsuperscript{-23} m\textsuperscript{2} at 10 MPa effective
pressure), indicating that the small amount of pre-existing pores are
predominantly unconnected (e.g., see Fig. 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
\selectlanguage{greek}μ\selectlanguage{english}m) and (ii) its virtually crack-free nature. The former is essential to
ensure a statistically significant number of grains (\textgreater{}10
grains per diameter) in the small (3 mm diameter x 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 micro-cracks, one sample was heated slowly to
600 \selectlanguage{ngerman}°C prior to deformation. Thermal stressing is one of the key
fracture-generating mechanisms in crustal rocks and is an effective
method for introducing micro-fractures into rock samples. Heating ACM to
elevated temperatures (\textgreater{}500 °C) induces significant,
permanent micro-crack damage, evident from photomicrographs (Meredith et
al., 2012) 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 micro-crack network with
average aperture of 0.3 \selectlanguage{greek}μ\selectlanguage{english}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 \selectlanguage{greek}μ\selectlanguage{english}CT data because this aperture
is less than one tenth the length of one pixel (2.7 \selectlanguage{greek}μ\selectlanguage{english}m). Due to the
partial volume effect, micro-cracks with an aperture smaller than half a
pixel are not visible (e.g., Voorn et al., 2013).

\subsection*{Experimental apparatus, sample assembly and loading
protocol}

{\label{experimental-apparatus-sample-assembly-and-loading-protocol}}

Synchrotron x-ray microtomography (\selectlanguage{greek}μ\selectlanguage{english}CT), in combination with x-ray
transparent pressure vessels (e.g., Fusseis et al., 2014; Renard et al.,
2016; Butler et al., 2017), allow the microstructural evolution of
deforming rock samples to be imaged directly, non-invasively
and\emph{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. \emph{In-situ} x-ray \selectlanguage{greek}μ\selectlanguage{english}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 \emph{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 (\textless{}1.3 kg) x-ray transparent tri-axial deformation
apparatus Mj\selectlanguage{ngerman}ö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 \selectlanguage{greek}μ\selectlanguage{english}CT rotation stage at the PSICHE beamline at SOLEIL
Synchrotron, Gif-sur-Yvette, France (Figure 2a,b). 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 \selectlanguage{greek}μ\selectlanguage{english}m. Even though this sample diameter is very small (required to
obtain high-resolution \selectlanguage{greek}μ\selectlanguage{english}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\textsuperscript{TM} high pressure
syringe pumps operated with QmixElements\textsuperscript{TM} software.

Experiments were conducted at room temperature under nominally dry
conditions. A reference \selectlanguage{greek}μ\selectlanguage{english}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 x
10\textsuperscript{-5} s\textsuperscript{-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 \selectlanguage{greek}μ\selectlanguage{english}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 \selectlanguage{greek}μ\selectlanguage{english}CT volume. Each full set of
scans was acquired in approx. 10 minutes. 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 \selectlanguage{greek}μ\selectlanguage{english}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, \(\sigma_{c}\)).

The differential stress is \(\sigma=\ \sigma_{1}-\ \sigma_{3}\), where \(\sigma_{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
\(\sigma_{3}\) is the radially-symmetric confining pressure. Axial
sample strain was calculated as\(\epsilon=\delta L/L_{0}\), where
\selectlanguage{greek}δ\selectlanguage{english}L is the change in length of the sample between the
starting \selectlanguage{greek}μ\selectlanguage{english}CT volume and the volume of interest and \(L_{0}\) is
the initial length of the sample. It was obtained directly from the
reconstructed \selectlanguage{greek}μ\selectlanguage{english}CT volumes by measuring the length change of the rock
core between two fixed locations in each volume.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image2/image2}

\end{center}
\end{figure}

\textbf{Figure 2:} Our x-ray transparent deformation apparatus, Mj\selectlanguage{ngerman}ölnir,
in (\textbf{a} ) schematic (Butler et al., 2017), and (\textbf{b} )
installed at the PSICHE beamline. (\textbf{c} ) Stepped loading
procedure. (\textbf{d} ) Reconstructed \selectlanguage{greek}μ\selectlanguage{english}CT image showing damage
accumulation at one end of the sample (this occurred in both samples) --
dashed white line shows the analyzed sub-volume.

\subsection*{X-ray imaging and image data
pre-processing}

{\label{x-ray-imaging-and-image-data-pre-processing}}

X-ray \selectlanguage{greek}μ\selectlanguage{english}CT data was acquired using an sCMOS (scientific
Complementary-Metal-Oxide Semiconductor) Hamamatsu ORCA Flash4.0 camera,
with a Rodagon 50 lens, giving about 2.5x magnification (effective pixel
size 2.7 \selectlanguage{greek}µ\selectlanguage{english}m), and a 250 \selectlanguage{greek}µ\selectlanguage{english}m thick LuAG:Ce scintillator. The white beam
with an average detected energy of about 66 keV was filtered with 1 mm
aluminium and 0.5 mm tungsten. During each scan, 1200 projections were
acquired over 180\selectlanguage{ngerman}°, 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 \selectlanguage{greek}µ\selectlanguage{english}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 3D volumes of 1700 x 1700 x
4102 equidimensional voxels, with a voxel edge length of 2.7 \selectlanguage{greek}μ\selectlanguage{english}m. These
volumes were then processed to extract the fracture network from the
reconstructed images. To deal efficiently with the huge size of each 3D
volume (approx. 40 GB) and speed up the subsequent processing, we
selected a sub-volume of interest -- the region in the failed samples
where the majority of damage had accumulated (Figure 2d, Table 1). Using
the Avizo\textsuperscript{TM} software package, this sub-volume was
extracted from each of the full 3D volumes and de-noised with an
anisotropic diffusion filter (stop value 0.4 over 4 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
datasets of manageable size (approx. 3 GB).

\textbf{Table 1:} Dimensions of the whole sample and analyzed
sub-volume, with uncertainties to two decimal places.\selectlanguage{english}
\begin{longtable}[]{@{}llllll@{}}
\toprule
\textbf{sample} & \textbf{sample diameter (mm)} & \textbf{sample length
(mm)} & \textbf{ROI length (mm)} & \textbf{ROI area
(mm\textsuperscript{2})} & \textbf{ROI volume
(mm\textsuperscript{3})}\tabularnewline
\midrule
\endhead
\textbf{untreated {[}ACfresh02{]}} & 2.96 & 8.73 & 3.33 & 6.62 &
22.03\tabularnewline
\textbf{heat-treated {[}ACHT01{]}} & 2.97 & 9.46 & 2.71 & 6.58 &
17.83\tabularnewline
\bottomrule
\end{longtable}

\subsection*{Segmentation}

{\label{segmentation}}

Each dataset of 3D grey-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 3D
image datasets. This is due to the range of greyscale values
accommodated by fractures of different apertures and the increasing
similarity of these grey-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 micro-cracks 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 3D volume (e.g., Descoteaux et
al\emph{.,} 2005). Attributes of this local curvature can be used to
distinguish planar features in the dataset (Text S1a in our Supporting
Information, SI). 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 3D
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).

\textbf{Table 2:} Input parameters for segmentation code. Definitions
given in Voorn et al. (2013).\selectlanguage{english}
\begin{longtable}[]{@{}llllllllll@{}}
\toprule
\textbf{smin} & \textbf{smax} & \textbf{sstep} & \textbf{blocksize} &
\textbf{avgmat} & \textbf{consthresh} & \textbf{maxmat} &
\textbf{padding} & \textbf{mincal} & \textbf{maxcal}\tabularnewline
\midrule
\endhead
2 & 8 & 1 & 1037 & 65535 & 0 & 56540 & 0 & 0.14118 &
0.32549\tabularnewline
\bottomrule
\end{longtable}

The binary segmented volumes were labelled and an intensity radius
threshold of 15 to 100 applied followed by a small spot filter to remove
round segmented objects with a radius of \textless{}50 pixels (134 \selectlanguage{greek}μ\selectlanguage{english}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 and Ruegsegger, 1997). Even with the segmentation
method described, there is still significant under-sampling 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).

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 (\(\sigma_{1}\)). The other two (\(x\) and
\(y\)) are the horizontal axes, which are perpendicular to
the loading direction and correspond to the confining stress
(\(\sigma_{2}=\sigma_{3}\)) with their directions arbitrarily assigned but
consistent between the two experiments. Void orientations are given in
terms of their dip \(\phi\) (deviation from horizontal) and
strike \(\theta\) (deviation from \(y\)).\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image3/image3}

\end{center}
\end{figure}

\par\null

\textbf{Figure 3:} (\textbf{a} ) Segmented fracture network (blue)
within a greyscale 3D \selectlanguage{greek}μ\selectlanguage{english}CT sub-volume at peak stress in the heat-treated
sample. (\textbf{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.

\subsection*{Analysis of the segmented microcrack
network}

{\label{analysis-of-the-segmented-microcrack-network}}

\subsubsection*{Porosity and the number of
fractures}

{\label{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
Avizo\textsuperscript{TM}. 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).

\subsubsection*{3D micro-crack orientations and
geometries}

{\label{d-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 Avizo\textsuperscript{TM}. 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).

\subsubsection*{Crack network scaling exponents and correlation
length}

{\label{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\selectlanguage{ngerman}é, 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 \selectlanguage{greek}µ\selectlanguage{english}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
Avizo\textsuperscript{TM}, 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 \selectlanguage{greek}µ\selectlanguage{english}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):

\begin{enumerate}
\tightlist
\item
  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.
\item
  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.
\item
  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\).
\end{enumerate}

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 \emph{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 r\textsuperscript{2} 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
r\textsuperscript{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 fractal dimension, \(D\), from
the relation\({P(L}_{i})\ \sim\ {L_{i}}^{D-1}\) \textbf{(} 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\) \selectlanguage{greek}μ\selectlanguage{english}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).

\section*{Results: Micro-fracture network
evolution}

{\label{results-micro-fracture-network-evolution}}

Here we present our 4D observations of the evolving segmented fracture
network in each sample, together with the distributions of micro-crack
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.

\subsection*{4D observations of micro-fractures and their
orientations}

{\label{d-observations-of-micro-fractures-and-their-orientations}}

The spatial evolution of micro-cracks differs significantly between the
two samples (Figure 4 and Figure 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
sub-parallel to \(\sigma_{1}\) (Figure 4F-O). Conversely, the
heat-treated sample failed along a localized shear fault, inclined at
30\selectlanguage{ngerman}° to \(\sigma_{1}\), along which pre-failure damage had already
accumulated (Figure 5K-P). In both samples, the localized damage zones
consisted mainly of shear and axial micro-fractures oriented between 5°
and 30° to \(\sigma_{1}\) (Figure 6a,b) 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\% strain and
64\% of the failure stress, \(\sigma_{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\%\(\sigma_{c}\). It is not clear exactly how the eventual fault
developed because the sample failed before the next scheduled image
time. The post-failure \selectlanguage{greek}µ\selectlanguage{english}CT sub-volume (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.2.1). However, close to peak stress (97\%
\(\sigma_{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\% \(\sigma_{c}\) (Table 4). Damage
progressively concentrated along this plane until localization around
the incipient fault plane became apparent (pink ellipse in Figure 5K) at
90\%\(\sigma_{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.2.1).

Microcrack dips, \(\phi\), became progressively more vertical
with increasing stress in both samples (Figure 6a,b), indicating the
preferential nucleation of tensile cracks. These cracks formed
en-echelon (Tapponnier and Brace, 1976; Kranz, 1979; Olson and Pollard,
1991; Reches and Lockner, 1994) and wing-crack \emph{(} Fairhurst and
Cook, 1966; Nemat-Nasser and Horii, 1982; Horii and Nemat-Nasser, 1985;
1986; Ashby and Hallam, 1986; Nemat-Nasser and Obata, 1988; Ashby and
Sammis, 1990) 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 onwards). In contrast, this process occurred
only in the heat-treated sample during localization around the eventual
fault plane (Figure 5L onwards), not during the initial localization
around the unfavorable conjugate. En-echelon and wing-crack arrays
formed at 1.45\(\pm\)0.01\% strain in both samples (Figure 4G
and Figure 5L). At this point the untreated sample was only at 70\%
\(\sigma_{c}\), compared to 90\% for the heat-treated sample.

One advantage of the 3D sampling enabled by \selectlanguage{greek}µ\selectlanguage{english}CT imaging is that we can
test the null hypothesis that the initial sample porosity is isotropic.
The optimal strikes of the segmented voids, \(\theta\), show a
predominant orientation in the initial porosity in both samples (Figure
6c,d). This starting anisotropy was more pronounced in the heat-treated
sample than in the untreated sample (33.0\selectlanguage{ngerman}±15.1\% compared with
14.3±11.8\% - see Table S1 in SI). Overall, anisotropy in the void
strike increased steadily throughout deformation in the heat-treated
sample but remained approximately constant in the untreated sample
(Table S1 in SI). 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 (\(\sigma_{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 in the SI), 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 micro-cracks
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\% \(\sigma_{c}\), bridging
four radial zones adjacent to the previous three, but in a conjugate
orientation. This bridging fault propagated up the sample (Video S2 in
the SI -- \(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 in the SI) shows that
localization initiated within the sample, not at the boundary, on the
site of a pre-existing void, precisely as anticipated by Griffith (1921,
1924). Subsequent micro-cracks 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 micro-cracks localized along it, before
propagating down the sample parallel to the \(z,x\) plane
(Video S4 in the SI -- \(y,z\)projection). Simultaneously,
micro-cracks 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, micro-cracks 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 micro-crack segments linked
by bridges of intact rock, along which the sample failed.

\textbf{Table 3:} Differential stress, \(\sigma\), and axial
strain, \(\epsilon\), for each \selectlanguage{greek}μ\selectlanguage{english}CT sub-volume in the untreated
sample {[}ACfresh02{]}. 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.\selectlanguage{english}
\begin{longtable}[]{@{}llllllllllllllllll@{}}
\toprule
\textbf{ACfresh02} & \textbf{A} & \textbf{B} & \textbf{C} & \textbf{D} &
\textbf{E} & \textbf{F} & \textbf{G} & \textbf{H} & \textbf{I} &
\textbf{J} & \textbf{K} & \textbf{L} & \textbf{M} & \textbf{N} &
\textbf{O} & \textbf{failure} & \textbf{P}\tabularnewline
\midrule
\endhead
\(\mathbf{\sigma}\) \textbf{(MPa)} & 7 & 41 & 60 & 79 & 97 & 116 & 126 &
135 & 144 & 154 & 158 & 163 & 168 & 172 & 177 & 182 & 127\tabularnewline
\(\mathbf{\epsilon}\) \textbf{(\%)} & 0.00 & 0.56 & 0.81 & 1.00 & 1.18 &
1.37 & 1.46 & 1.62 & 1.71 & 1.77 & 1.84 & 1.87 & 1.93 & 1.99 & 2.08 & -
& 25.79\tabularnewline
\bottomrule
\end{longtable}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image4/image4}

\end{center}
\end{figure}

\textbf{Figure 4:} Projection in the \(y,z\) plane of
segmented microfractures for the untreated sample {[}ACfresh02{]}.
Orange ellipses show the onset of localization, pink ellipses indicate
formation of additional localized zone, completing the radial pattern.
Letters refer to the stress-strain steps in Table 3. Color-coding shows
the local fracture aperture at each voxel. See Figure S5 in SI for zip
file of high resolution images.

\textbf{Table 4:} Differential stress, \(\sigma\), and axial
strain, \(\epsilon\), for each \selectlanguage{greek}μ\selectlanguage{english}CT sub-volume in the heat-treated
sample {[}ACHT01{]}. 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.\selectlanguage{english}
\begin{longtable}[]{@{}lllllllllllllllll@{}}
\toprule
\textbf{ACHT01} & \textbf{A} & \textbf{B} & \textbf{C} & \textbf{D} &
\textbf{E} & \textbf{F} & \textbf{G} & \textbf{H} & \textbf{I} &
\textbf{J} & \textbf{K} & \textbf{L} & \textbf{M} & \textbf{N} &
\textbf{O (failure)} & \textbf{P}\tabularnewline
\midrule
\endhead
\(\mathbf{\sigma}\) \textbf{(MPa)} & 7 & 41 & 59 & 78 & 97 & 115 & 125 &
134 & 153 & 162 & 167 & 171 & 176 & 181 & 185 & 102\tabularnewline
\(\mathbf{\epsilon}\) \textbf{(\%)} & 0.03 & 0.35 & 0.63 & 0.84 & 0.92 &
1.04 & 1.09 & 1.24 & 1.30 & 1.38 & 1.41 & 1.44 & 1.53 & 1.73 & 1.87 &
5.65\tabularnewline
\bottomrule
\end{longtable}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image5/image5}

\end{center}
\end{figure}

\textbf{Figure 5:} Projection in the \(y,z\) plane of
segmented microcracks for the heat-treated sample {[}ACHT01{]}. Orange
ellipse indicates the onset of damage localization, pink ellipse
indicates localization along the emerging fault plane. Letters refer to
the stress-strain steps in Table 4. Color-coding shows the local
fracture aperture at each voxel. See Figure S6 in SI for zip file of
high resolution images.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image6/image6}

\end{center}
\end{figure}

\textbf{Figure 6:} Void (\textbf{a,b} ) dip \(\phi\) and
(\textbf{c,d} ) strike \(\theta\)probability distributions at
each strain increment for (\textbf{a,c} ) the untreated sample and
(\textbf{b,d} ) the heat-treated sample, calculated using the
object-based best-fitting ellipsoid approach. Inset shows the measured
angles with respect to the principal stress axes:\(\phi=90\) is
parallel to the \(z\)-axis, while \(\theta=0\) is
parallel with the \(y\)-axis. Labels on the dip
distributions show the angles at which different damage mechanisms are
prevalent. Strike distributions are normalized to the minimum
\(P\left(\theta\right)\)of the starting porosity. Blue distributions
overlaying the rest are post-failure.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image7/image7}

\end{center}
\end{figure}

\textbf{Figure 7:} Slice through the sample in the \(y,z\)
projection showing in (\textbf{a} ) the untreated sample close to peak
stress (Figure 4O) and (\textbf{b} ) the heat-treated sample at peak
stress (Figure 5O) the formation of tensile (\textbf{c} ) en-echelon
(Olson and Pollard, 1991; Thomas and Pollard, 1993) and (\textbf{d} )
wing-cracks (Fairhurst and Cook, 1966; Fig. 15 in Damjanac and
Fairhurst, 2010).\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image8/image8}

\end{center}
\end{figure}

\textbf{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/\sqrt{\pi(9+n)}\), 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 micro-crack 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 micro-cracks
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).

\subsection*{Evolution of specimen and microcrack characteristics with
strain}

{\label{evolution-of-specimen-and-microcrack-characteristics-with-strain}}

\subsubsection*{Stress, porosity and the number of
voids}

{\label{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\% \(\sigma_{c}\) in the heat-treated sample (Figure 5M), but
only from 97\% \(\sigma_{c}\) in the untreated sample (Figure 4O).

The heat-treated sample had lower pre-existing porosity than the
untreated sample (\(\varphi_{0\ HT}=0.62\varphi_{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 (Figure 9b,c).
However, this observation only accounts for voids visible above the
detection threshold of the segmentation algorithm (a void volume of
3000-4000 \selectlanguage{greek}μ\selectlanguage{english}m\textsuperscript{3} -- see Section 3.3.1), and does not
include unresolved nano-scale 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 ten-fold increase in porosity, \(\varphi,\)
over the duration of their respective deformation experiment (Figure
9b), but only a two-fold increase in the total number of voids,
\(N\), in the heat-treated sample, compared with a nearly
three-fold 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 (Figure 4E-F) is evident as a
large jump in \(N\) of 600 voids at 1.37\% strain, and a
corresponding three-fold increase in\(\varphi\) (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 \(\varphi_{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 \(\varphi\) 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 towards failure in both \(\varphi\) 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 micro-cracks 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 \(\varphi\) and
\(N\) with strain is best described with simple power-law
models (Figure 9b,c); i.e., they have the lowest AICc,
(Tables S2 and S3 in SI). The exponent is 3.1 for both variables in the
untreated sample, compared with 8.8 and 7.7 for \(\varphi\) and
\(N\) respectively in the heat-treated sample, showing an
acceleration towards 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 \(\varphi\) 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\(\varphi\), confirming that crack growth played
an increasingly important role closer to failure.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image9/image9}

\end{center}
\end{figure}

\textbf{Figure 9:} Evolution of (\textbf{a} ) differential stress,
\(\sigma\), (\textbf{b} ) porosity,\(\varphi\) and
(\textbf{c} ) the number of segmented voids, \(N\), with
increasing axial strain, \(\epsilon\) 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 \selectlanguage{greek}μ\selectlanguage{english}CT volumes (Figure 4
and Figure 5). Solid lines in (\textbf{a} ) show the region of data used
to calculate the Young's moduli (11.3 and 9.5 GPa respectively). Solid
lines in (\textbf{b} ) and (\textbf{c} ) show the preferred simple
power-law models with 95\% confidence intervals shown as dashed lines;
see text for exponents.

\subsubsection*{Micro-crack geometry}

{\label{micro-crack-geometry}}

To establish empirically how the micro-crack 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 \selectlanguage{greek}µ\selectlanguage{english}CT sub-volume (Figure 10a,b). We also show the mean
ellipsoid eigenvalue ratios, used to infer the evolution of void aspect
ratio (Figure 10c,d). We present two ratios: (i) the smallest to the
medium eigenvalue of the covariance matrix (Section 2.5.2 and Text S1b
in SI), 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, \(\overset{\overline{}}{r}\), were about the same
size in both samples. In the untreated sample (Figure
10a)\({\overset{\overline{}}{r}}_{\text{major}}\) and\({\overset{\overline{}}{r}}_{\text{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
micro-cracks localized along the optimally oriented damage zone (Figure
5K), after the onset of localization. In both samples,\(r_{\text{major}}\)
were oriented approximately parallel to the strike of the eventual fault
plane (Figure 6c,d), with their mean values,\({\overset{\overline{}}{r}}_{\text{major}}\),
increasing more quickly than \({\overset{\overline{}}{r}}_{\text{medium}}\), showing that
micro-cracks grew twice as fast along strike (i.e., perpendicular
to\(\sigma_{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
(\({\overset{\overline{}}{r}}_{\text{medium}}\)) increased from 2.5 to 4 times the crack aperture
(\({\overset{\overline{}}{r}}_{\text{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 micro-cracks 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.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image10/image10}

\end{center}
\end{figure}

\textbf{Figure 10:} (\textbf{a,b} ) Evolution with strain of the mean
void ellipsoid radii, \({\overset{\overline{}}{r}}_{\text{major}}\) (blue
circles),\({\overset{\overline{}}{r}}_{\text{medium}}\) (orange circles) and\({\overset{\overline{}}{r}}_{\text{minor}}\)
(yellow circles), for the (\textbf{a} ) untreated and (\textbf{b} )
heat-treated samples. (\textbf{c,d} ) Evolution of mean void eigenvalue
ratios (\textbf{c} )\({\overset{\overline{}}{r}}_{\text{med}}/{\overset{\overline{}}{r}}_{\text{maj}}\)and (\textbf{d}
)\({\overset{\overline{}}{r}}_{\min}/{\overset{\overline{}}{r}}_{\text{med}}\)with strain in the untreated (blue circles) and
heat-treated (orange circles) samples. Voids become flatter or more
elongate as the respective ratio \(\rightarrow 0\). Error bars show the
standard error of the mean in each \selectlanguage{greek}μ\selectlanguage{english}CT sub-volume. Dash-dot lines show
the failure strain for each sample while dashed lines show the onset of
localization as seen in the \selectlanguage{greek}μ\selectlanguage{english}CT volumes.

\subsection*{Evidence for phase transition
style}

{\label{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
\(\sigma\)and axial sample strain \(\epsilon\). 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, \(\xi\)(linear dimension of the largest void) evolves
as a function of stress. We then analyze the evolution of
\(\xi\), \(\beta\) and \(D\) (the
void volume and inter-void length exponents respectively, defined in
Section 2.5.3) as a function of strain.

\subsubsection*{Microcrack volume and inter-crack length
distributions}

{\label{microcrack-volume-and-inter-crack-length-distributions}}

Both samples show an approximately power-law complementary probability
distribution in void volume, \(V_{i}\), (Figure 11a,b), 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}\) (Figure
11c,d), within a finite range, identified as \(30<L_{i}<1350\) \selectlanguage{greek}μ\selectlanguage{english}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 \(\beta\) 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 3000 to 4000 \selectlanguage{greek}μ\selectlanguage{english}m\textsuperscript{3}, roughly
equivalent to a void aperture of 14-16 \selectlanguage{greek}μ\selectlanguage{english}m. This is much larger than the
theoretical detection threshold of half the pixel size (1.3 \selectlanguage{greek}μ\selectlanguage{english}m)
consistent with under-sampling of very narrow cracks during
segmentation. Void volumes in the untreated sample are best described
(i.e., have the lowestBIC) 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\% \(\sigma_{c}\), two stages
before the onset of localization (Figure S2 in the SI). In contrast,
void volumes in the heat-treated sample are best described by the GR
distribution throughout the experiment.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image11/image11}

\end{center}
\end{figure}

\textbf{Figure 11:} (\textbf{a,b} ) Cumulative complementary probability
distributions of void volumes, \(V_{i}\), computed with a bin
width of 4000 \selectlanguage{greek}μ\selectlanguage{english}m\textsuperscript{3}, and (\textbf{c,d} ) cumulative
probability distributions of inter-void lengths, \(L_{i}\),
computed with 256 bins, at each increment of stress and strain in the
(\textbf{a,c} ) untreated and (\textbf{b,d} ) heat-treated samples. As
described in Section 2.5.3, we obtained maximum likelihood estimates for
\(\beta\) from the cumulative complementary \(V_{i}\)
data (see Figure S1 in SI) using the models of Kagan (2002), and fitted
the incremental \(L_{i}\) data using linear regression in
log-log space (see Figure S3 in SI) to find \(D\), after
Turcotte (1997). Blue distributions overlaying the rest are
post-failure.

\subsubsection*{An inverse power-law acceleration to
failure?}

{\label{an-inverse-power-law-acceleration-to-failure}}

Parameters for an inverse power-law acceleration to failure for the
normalized correlation length, \(\xi/l\), were obtained for
both samples (Figure 12), using data observed in segmented \selectlanguage{greek}μ\selectlanguage{english}CT volumes
between Figure 4F-O and Figure 5H-O 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,
\(\sigma_{p}\), is far from the observed failure stress,
\(\sigma_{c}\). The likelihood that the inverse power-law model fits
the data as well as the exponential model is just 3\% (Table S5 in SI).
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,\(\sigma_{p}\) is
accurate to within 3\% of \(\sigma_{c}\), while the asymptote of
the exponential model is further from \(\sigma_{c}\).The likelihood
of the exponential model fitting the data as well as the inverse
power-law is 72\% (Table S5 in SI), which although relatively high, is
not significant (\textgreater{}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.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image12/image12}

\end{center}
\end{figure}

\textbf{Figure 12:} Evidence for an inverse power-law acceleration
(solid lines) with respect to stress in the correlation length
(normalized by the length of the analyzed sub-volume):
\(\xi/l=k{(\sigma_{p}-\sigma)}^{-p}\) in the heat-treated sample (orange), to a predicted
failure stress,\(\sigma_{p}\), within 3\% of the observed failure
stress,\(\sigma_{c}\). This compares with the same model yielding a
much poorer prediction in the untreated sample (blue),
where\(\sigma_{p}=2.4\ \sigma_{c}\). Dotted lines correspond to the exponential
model, \(\xi/l=h\ \exp(q\sigma)\), 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 in SI, along with the
AICc values). The dash-dot lines show
\(\sigma_{c}\)for both samples and the dashed line shows
\(\sigma_{p}\) for the heat-treated sample.

\subsubsection*{Evolution of crack population metrics with respect to
strain}

{\label{evolution-of-crack-population-metrics-with-respect-to-strain}}

Both samples show a systematic increase in the correlation
length\(\xi\) towards failure as a function of strain, with
failure occurring when \(\xi\) increased beyond 200 \selectlanguage{greek}µ\selectlanguage{english}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 \selectlanguage{greek}μ\selectlanguage{english}m -- see Section 2.1). This implies that the sample
breaks when whole grains break. The nature of the increase in
\(\xi\left(\epsilon\right)\) in the untreated sample (blue) was exponential
(i.e., had the lowestAICc -- see Table S4 in SI) up to a
finite \(\xi\) that fluctuated around 200 \selectlanguage{greek}µ\selectlanguage{english}m before failure.
Conversely, in the heat-treated sample (orange) \(\xi\left(\epsilon\right)\)
preferred a simple power-law acceleration (Figure S4 and Table S4 in SI)
to failure, failing when \(\xi>\) 200 \selectlanguage{greek}μ\selectlanguage{english}m. The exponent (6.9)
is the same (within error) as the exponents for the evolution of
\(\varphi\) and \(N\) with strain (Section 3.2.1;
Tables S2-S4 in the SI), independently confirming that in this sample
crack growth played an increasingly important role closer to failure.
The power-law acceleration emerges only once fractures began to localize
along the optimally oriented damage zone (at 90\% \(\sigma_{c}\)
and \(\xi>\ \xi_{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 \(\beta\) for
their volume distributions (Figure 13b). In the untreated sample
(blue)\(\beta\) 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
(Figure 9b,c). This trade-off is consistent with independent
observations from acoustic emissions (Graham et al., 2010) and models
(e.g. Brantut et al., 2012; 2014) of the competition between compaction
and dilatancy during the quasi-elastic region of the stress-strain curve
(Figure 9a). Beyond this peak, \(\beta\) 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 (Figure 9c and Figure 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) \(\beta\)
decreased throughout, reflecting an increase in the relative proportion
of larger micro-cracks. 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 \(\beta\) is a
clear precursor to failure, corresponding to propagation of the shear
zone through the sample (Figure 5M-O). 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 in the two time windows
before the onset of localization, 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 micro-cracks 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
micro-crack 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.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image13/image13}

\end{center}
\end{figure}

\textbf{Figure 13:} Evolution with strain for (\textbf{a} ) correlation
length \(\xi\), (\textbf{b} ) microcrack volume exponent
\(\beta\), and (\textbf{c} ) fractal
dimension\(D\), for the untreated (blue) and heat-treated
(orange) samples. Error bars show (\textbf{b} ) the standard error
(\(\beta/\sqrt{N}\)) of the maximum likelihood estimate of
\(\beta\), where \(N\) is the number of cracks,
and (\textbf{c} ) the 95\% confidence intervals of the linear regression
fit for \(D\). Dash-dot lines show the strain at which each
sample failed, while dashed lines show the onset of localization as seen
in the \selectlanguage{greek}μ\selectlanguage{english}CT volumes.

\section*{Discussion}

{\label{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 (\selectlanguage{greek}μ\selectlanguage{english}CT) observations
of \emph{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.

\subsection*{Microcrack network
evolution}

{\label{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 \emph{in-situ} \selectlanguage{greek}µ\selectlanguage{english}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 \textgreater{}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
recent\emph{in-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 \selectlanguage{ngerman}±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.

\subsection*{Scaling, phase transition style and predictability of
failure
time}

{\label{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 \selectlanguage{greek}µ\selectlanguage{english}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\selectlanguage{ngerman}è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
stress\emph{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 and Hauptfleisch, 1999), with a homogeneous
microstructure (Oesterling, 2004), a very low permeability
(10\textsuperscript{-19}m\textsuperscript{2}) 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 \selectlanguage{greek}μ\selectlanguage{english}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 \selectlanguage{greek}μ\selectlanguage{english}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 \selectlanguage{greek}µ\selectlanguage{english}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.

\subsection*{Suggestions for future
work}

{\label{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.

\section*{Conclusions}

{\label{conclusions}}

Our \emph{in-situ} time-resolved x-ray \selectlanguage{greek}µ\selectlanguage{english}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 micro-cracks have power-law
frequency-volume and inter-crack length distributions over a finite
scaling range, irrespective of the degree of starting heterogeneity,
with well-determined scaling exponents \(\beta\) (the size
exponent) and\(D\) (the correlation dimension). The
inferred correlation length\(\xi\) 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 \(\beta\) 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.

\section*{Acknowledgments, Samples, and
Data}

{\label{acknowledgments-samples-and-data}}

This work is supported by the UK's Natural Environment Research Council
(NERC) through the CATFAIL project NE/R001693/1 `Catastrophic failure:
what controls precursory localization in rocks?' We acknowledge the
beamline PSICHE at SOLEIL for provision of synchrotron radiation
facilities (standard proposal 20160434) and thank our reviewers Jess
McBeck and Phil Benson for their thoughtful reviews and suggestions. We
would also like to thank the University of Edinburgh Geosciences
Workshop for their support in developing the experimental apparatus, and
Andy Bell for useful discussions about modelling the correlation length
evolution. We declare no conflict of interest. The data supporting our
conclusions can be found in main text, the Supporting Information and in
the \selectlanguage{greek}μ\selectlanguage{english}CT datasets held at the NERC repository
(http://data.ceda.ac.uk/ngdc/R001693-1). When using the \selectlanguage{greek}μ\selectlanguage{english}CT datasets,
please cite Cartwright-Taylor et al. (2020) and see
https://doi.org/10.5285/0dc00069-8da8-474a-8993-b63ef5c25fb8 for the
metadata.

\section*{Appendix A. Supporting Information
(SI)}

{\label{appendix-a.-supporting-information-si}}

Supporting information related to this article can be found in the
accompanying document.

\section*{References}

{\label{references}}

Abercrombie, R.E. (1995), Earthquake source scaling relationships from
-1 to 5 ML using seismograms recorded at 2.5-km depth, J. Geophys. Res.
100, 24015--24036.

Aki, K. (1965), Maximum likelihood estimate of \(b\) in the
formula\(\log N=a-bM\) and its confidence limits, Bull. Earthquake
Res. Inst. 43, 237-239.

Alava, M.J., Nukala, P.K.V. and Zapperi, S. (2008), Fracture size
effects from disordered lattice models, Int. J. Fract. 154, 51-59.

Alber, M. and Hauptfleisch, U. (1999), Generation and visualization of
micro-fractures in Carrara marble for estimating fracture toughness,
fracture shear and fracture normal stiffness, Int. J. Rock Mech. Min.
Sci. 36, 1065-1071.

Andrew, M. (2018), A quantified study of segmentation techniques on
synthetic geological XRM and FIB-SEM images, Computational Geosciences
22, 1503-1512.

Ashby M.F. and Sammis, C.G. (1990), The damage mechanics of brittle
solids in compression, PAGEOPH 133, 489-521.

Ashby, M.F. and Hallam, S.D. (1986), The failure of brittle solids
containing small cracks under compressive stress states, Acta Metall.
34, 497-510.

Bandini, A., Berry, P., Bemporad, E., and Sebastiani, M. (2012), Effects
of intra-crystalline microcracks on the mechanical behavior of a marble
under indentation, Int. J. Rock Mech. Min. Sci. 54, 47-55.

Bell, A.F., Naylor, M. and Main, I.G. (2013a), Convergence of the
frequency-size distribution of global earthquakes, Geophys. Res. Lett.
40, 2585-2589.

Bell, A.F., Naylor, M. and Main, I.G. (2013b), The limits of
predictability of volcanic eruptions from accelerating rates of
earthquakes, Geophys. J. Int. 194, 1541-1553.

Bonnet, E., Bour, O., Odling, N.E., Davy, P., Main, I., Cowie, P. and
Berkowitz, B. (2001), Scaling of fracture systems in geological media,
Rev. Geophys. 39, 347-383.

Bouchaud, E. (1997), Scaling properties of cracks, J. Phys. Condens.
Matter 9, 4319-4344.

Bouchaud, E., Lapasset, G., and Plan\selectlanguage{ngerman}ès, J. (1990), Fractal dimension of
fractured surfaces: a universal value?, Europhys. Lett. 13, 73-79.

Brace, W.F., Paulding, B. and Scholz, C. (1966), Dilatancy in the
fracture of crystalline rocks. J. Geophys. Res. 71, 3939-3953.

Brantut, N., Baud, P., Heap, M.J. and Meredith, P.G. (2012),
Micromechanics of brittle creep in rocks, J. Geophys. Res. 117, B08412.

Brantut, N., Heap, M.J., Baud, P. and Meredith, P.G. (2014), Rate- and
strain-dependent brittle deformation of rocks, J. Geophys. Res. Solid
Earth 119, 1818-1836.

Bruce, A. and D. Wallace (1989). Critical point phenomena: universal
physics at large length scales, in: Davies, P. (ed.), The New Physics,
Cambridge University Press, Cambridge UK.

Bruner, W.M. (1984), Crack growth during unroofing of crustal rocks:
effects of thermoelastic behavior and near-surface stresses, J. Geophys.
Res. 89, 4167-4184.

Bruner, W.M. (1979), Crack growth and the thermoelastic behavior of
rocks, J. Geophys. Res. 84, 5578-5590.

Bufe, C.G. and Varnes, D.J. (1993), Predictive modeling of the seismic
cycle of the greater San Francisco Bay region, J. Geophys. Res. 98,
9871-9883.

Burnham, K.P. and Anderson, D.R. (2002), Model selection and multi-model
inference: a practical information-theoretic approach (2nd ed.),
Springer, New York.

Butler, I.B., Flynn, M., Fusseis, F. and Cartwright-Taylor, A. (2017),
Mjolnir: a novel x-ray transparent triaxial rock deformation apparatus,
ICTMS2017-56, 3rd International Conference on Tomography of Materials
and Structures, Lund, Sweden, 26-30 June.

Cartwright-Taylor, A. (2015), Deformation-induced electric currents in
marble under simulated crustal conditions: non-extensivity,
superstatistical dynamics and implications for earthquake hazard, PhD
thesis, University College London.

Cartwright-Taylor, A., Main, I.G., Butler, I.B., Fusseis, F., Flynn M.
and King, A. (2020): In-situ rock deformation and micron-scale crack
network evolution: a high-resolution time-resolved x-ray
micro-tomography dataset, British Geological Survey (Dataset)
https://doi.org/10.5285/0dc00069-8da8-474a-8993-b63ef5c25fb8.

Ceva, H. and Perazzo, R.P.J. (1993), From self-organized criticality to
first-order-like behaviour: A new type of percolative transition, Phys.
Rev. E 48, 157-160.

Cho, N., Martin, C.D. and Sego, D.S. (2007), A clumped particle model
for rock, Int. J. Rock Mech. Min. Sci. 44, 997-1010.

Clint, O.C., Meredith, P.G. and Main, I.G. (2001), Relation between
crack damage and permeability near the percolation threshold in a near
perfect crystalline rock, Geophys. Res. Abstr. 3, 346.

Cornelius, R.R. and Voight, B. (1994), Seismological aspects of the
1989-1990 eruption at Redoubt Volcano, Alaska: The Materials Failure
Forecast Method (FFM) with RSAM and SSAM seismic data, J. Volcanol.
Geotherm. Res. 62, 469-498.

Cowie, P.A., Sornette, D. and Vanneste, C. (1995), Multifractal scaling
properties of a growing fault population, Geophys. J. Int. 122, 457-469.

Crawford, B.R., Smart, B.G.D., Main, I.G. and Liakopoulou-Morris, F.
(1995). Strength characteristics and shear acoustic anisotropy of rock
core subjected to true triaxial compression, Int. J. Rock Mech. Min.
Sci. 32, 189-200.

Damjanac, B. and Fairhurst, C. (2010), Evidence for a long-term strength
threshold in crystalline rock, Rock Mech. Rock Eng. 43, 513-531.

David, C., Menendez, B. and Darot, M. (1999), Influence of
stress-induced and thermal cracking on physical properties and
microstructure of La Peyratte granite, Int. J. Rock Mech. Min. Sci. 36,
433--448.

Descoteaux, M., Audette, M., Chinzei, K. and Siddiqi,K. (2005), Bone
enhancement filtering: application to sinus bone segmentation and
simulation of pituitary surgery, in: Duncan, J.S. and Gerig, G. (eds.),
Proceedings of Medical Image Computing and Computer-Assisted
Intervention (MICCAI), Palm Springs, California, USA, Lecture Notes in
Computer Science 3749, 09--16.

Desrues, J., Chambon, R., Mokni, M. and Mazerolle, F. (1996), Void ratio
evolution inside shear bands in triaxial sand specimens studied by
computed tomography, Geotechnique 46, 529-546.

Fairhurst, C. and Cook, N.G.W. (1966), The phenomenon of rock splitting
parallel to the direction of maximum compression in the neighborhood of
a surface, In: Zeitlen, J.G. (ed) Proceedings of First Congress of the
International Society of Rock Mechanics (Lisbon, September-October 1966)
Vol 1., LNEC, Lisbon, 687-692.

Faulkner, D.R., Mitchell, T.M., Healy, D. and Heap, M.J. (2006), Slip on
`weak' faults by the rotation of regional stress in the damage zone,
Nature 444, 922-925.

Fredrich, J.T., and Wong T.-F. (1986), Micromechanics of thermally
induced cracking in three crustal rocks, J. Geophys. Res. 91,
12743--12764.

Fredrich, J.T., Evans, B., and Wong, T.-F. (1989), Micromechanics of the
brittle to plastic transition in Carrara marble. J. Geophys. Res. Solid
Earth 94, 4129-4145.

Fusseis, F., Steeb, H., Xiao, X., Zhu, W., Butler, I., Elphick, S. and
Mäder, U. (2014), A low-cost X-ray-transparent experimental cell for
synchrotron-based X-ray microtomography studies under geological
reservoir conditions, J. Synchrotron Rad. 21, 251-253.

Girard, L., Amitrano, D. and Weiss, J. (2010), Failure as a critical
phenomenon in a progressive damage model, J. Stat. Mech. Theory Exp.
P01013.

Goodfellow, S.D. and Young, R.P. (2014), A laboratory acoustic emission
experiment under in-situ conditions, Geophys. Res. Lett. 41, 3422-3430.

Graham, C.C., Stanchits, S., Main, I.G. and Dresen, G. (2010),
Comparison of polarity and moment tensor inversion methods for source
analysis of acoustic emission data, Int. J. Rock Mech. Min. Sci. 47,
161-169.

Greenhough, J. and Main, I.G. (2008), A Poisson model for earthquake
frequency uncertainties in seismic hazard analysis, Geophys. Res. Lett.
35, L19313.

Griffith, A.A. (1921), The phenomena of rupture and flow in solids,
Phil. Trans. R. Soc. A 221, 163-198.

Griffith, A.A. (1924), Theory of rupture, Proc. First Int. Congress
Appl. Mech., Delft, 55-63.

Griffiths, L., Lengliné, O., Heap, M.J., Baud, P. and Schmittbuhl, J.
(2018), Thermal cracking in Westerly granite monitored using direct wave
velocity, coda wave interferometry and acoustic emissions, J, Geophys.
Res. Solid Earth 123, 2246-2261.

Guéguen, Y. and Schubnel, A. (2003), Elastic wave velocities and
permeability of cracked rocks, Tectonophys. 370, 163-176.

Hallbauer, D.K., Wagner, H. and Cook N.G.W. (1973), Some observations
concerning the microscopic and mechanical behaviour of quartzite
specimens in stiff, triaxial compression tests, Int. J. Rock Mech. Min.
Sci. Geomech. Abstr. 10, 713-726.

Haneberg, W.C. (2004), Computational Geosciences with Mathematica,
Springer-Verlag, Berlin Heidelberg New York.

Healy, D., Blenkinsop, T.G., Timms, N.E., Meredith, P.G., Mitchell, T.M.
and Cooke, M.L. (2015), Polymodal faulting: time for a new angle on
shear failure, J. Struct. Geol. 80, 57-71.

Hildebrand, T. and Rüegsegger, P. (1997), A new method for the
model\selectlanguage{english}-independent assessment of thickness in three\selectlanguage{english}-dimensional images,
J. Microscopy 185, 67-75.

Horii, H. and Nemat-Nasser, S. (1985), Compression-induced microcrack
growth in brittle solids: axial splitting and shear failure, J. Geophys.
Res. 90, 3105-3125.

Horii, H. and Nemat-Nasser, S. (1986), Brittle failure in compression:
splitting, faulting and brittle-ductile transition, Phil. Trans. R. Soc.
A 319, 337-374.

Hurvich, C.M. and Tsai, C.-L. (1989), Regression and time series model
selection in small samples, Biometrika 76, 297-307.

Kagan, Y.Y. (2002), Seismic moment distribution revisited: i.
statistical results, Geophys. J. Int. 148, 520-541.

Kagan, Y.Y. and Schoenberg, F. (2001), Estimation of the upper cutoff
parameter for the tapered Pareto distribution, J. Appl. Probab. 38,
168-185.

Kamb, W.B. (1959), Ice petrofabric observations from Blue Glacier,
Washington, in relation to theory and experiment, J. Geophys. Res. 64,
1891-1909.

Kandula, N., Cordonnier, B., Boller, E., Weiss, J., Dysthe, D.K. and
Renard, F. (2019), Dynamics of microscale precursors during brittle
compressive failure in Carrara marble, J. Geophys. Res. Solid Earth,
124, 6121-6139.

Kijko, A. and Graham, G. (1998), Parametric-historic procedure for
probabilistic seismic hazard analysis part i: estimation of maximum
regional magnitude \(m_{\max}\), PAGEOPH 152, 413-442.

Kilburn, C.R.J. (2003), Multiscale fracturing as a key to forecasting
volcanic eruptions, J. Volcanol. Geotherm. Res. 125, 271-289.

Kilburn, C.R.J. and Voight, B. (1998), Slow rock fracture as eruption
precursor at Soufriere Hills, Geophys. Res. Lett. 25, 3665-3668.

Kosterlitz, J.M. (1974), The critical properties of the two-dimensional
xy model, J. Phys. C: Solid State Phys. 7, 1046-1060.

Kosterlitz, J.M. and Thouless, D.J. (1973), Ordering, metastability and
phase transitions in two-dimensional systems, J. Phys. C: Solid State
Phys. 6, 1181-1203.

Kranz, R.L. (1979), Crack-crack and crack-pore interactions in stressed
granite, Int. J. Rock Mech. Min. Sci. Geomech. Abstr. 16, 37-47.

Kun, F., Pal, G., Varga, I. and Main, I.G. (2018), Effect of disorder on
the spatial structure of damage in slowly compressed porous rocks, Phil.
Trans. R. Soc. A 377, 20170393.

Kun, F., Varga, I., Lennartz-Sassinek, S. and Main, I.G. (2013),
Approach to failure in porous granular materials under compression,
Phys. Rev. E 88, 062207.

Lajtai, E.Z. (1974), Brittle fracture in compression. Int. J. Fract. 10,
525-536.

Lei, X. and Satoh, T. (2007), Indicators of critical point behavior
prior to rock failure inferred from pre-failure damage, Tectonophys.
431, 97-111.

Lei, X., Kusunose, K., Nishizawa, O., Cho, A. and Satoh, T. (2000), On
the spatiotemporal distribution of acoustic emissions in two granitic
rocks under triaxial compression: the role of pre-existing cracks,
Geophys. Res. Lett. 27, 1997-2000.

Lennartz-Sassinek, S., Main, I.G., Zaiser, M. and Graham, C.C. (2014),
Acceleration and localization of subcritical crack growth in a natural
composite material, Phys. Rev. E 90, 052401.

Leonard, T. and Hsu, J.S.J. (1999), Bayesian Methods, Cambridge
University Press, Cambridge UK.

Liakopoulou-Morris, F., Main, I.G., Crawford, B.R. and Smart, B.G.D.
(1994), Microseismic properties of a homogeneous sandstone during fault
nucleation and frictional sliding, Geophys. J. Int. 119, 219-230.

Lockner, D., Byerlee, J.D., Kuksenko, V., Ponomarev, A. and Sidorin, A.
(1992), Observations of quasi-static fault growth from acoustic
emissions, in: Evans, B. and Wong, T.-F. (eds.), Fault Mechanics and
Transport Properties of Rocks (1st ed.), International Geophysics 51,
3-31.

Lockner, D., Byerlee, J.D., Kuksenko, V., Ponomarev, A. and Sidorin, A.
(1991), Quasi-static fault growth and shear fracture energy in granite.
Nature 350, 39-42.

Lominitz-Adler, J., Knopoff, L. and Martinez-Mekler, G. (1992),
Avalanches and epidemic models of fracturing in earthquakes, Phys. Rev.
A 45, 2211-2221.

Main, I.G. (1996), Statistical physics, seismogenesis and seismic
hazard, Rev. Geophys. 34, 433-462.

Main, I.G., Meredith, P.G., Sammonds, P.R. and Jones, C. (1990),
Influence of fractal flaw distributions on rock deformation in the
brittle field, in: Knipe, R.J. and Rutter, E.H. (eds.), Deformation
Mechanisms, Rheology and Tectonics, Geol. Soc. Spec. Pub. 54, 81-96.

Main, I.G. (1992). Damage mechanics with long-range interactions:
correlation between the seismic b-value and the two point correlation
dimension, Geophys. J. Int. 111, 531-541.

Main, I.G., Sammonds, P.R. and Meredith, P.G. (1993), Application of a
modified Griffith criterion to the evolution of fractal damage during
compressional rock failure, Geophys. J. Int. 115, 367-380.

Mandelbrot, B., Passoja, D.E., and Paullay, A.J. (1984), Fractal
character of fracture surfaces of metals, Nature 308, 721-722.

Meijering, E.H.W. (2010), FeatureJ 1.6.0, Biomedical Imaging Group
Rotterdam, Erasmus MC, University Medical Center Rotterdam, The
Netherlands, 2002-2010
(http://www.imagescience.org/meijering/software/featurej/).

Meredith, P.G., Clint, O.C., Ngwenya, B., Main, I.G., Odling, N.W.A. and
Elphick, S.C., (2005), Crack damage and permeability evolution near the
percolation threshold in a near-perfect crystalline rock, in: Shaw,
R.P., (ed.) Understanding the Micro to Macro Behavior of Rock-Fluid
Systems, Geol. Soc. Spec. Pub. 249, 159-160.

Meredith, P.G., Main, I.G., Clint, O.C. and Li, L. (2012), On the
threshold of flow in a tight natural rock, Geophys. Res. Lett. 39,
L04307.

Mignan, A. and Woessner, J. (2012), Estimating the magnitude of
completeness for earthquake catalogues. Community Online Resource for
Statistical Seismicity Analysis, pp. 1--45.

Mirone, A., Brun, E., Gouillart, E., Tafforeau, P. and Kieffer, J.
(2014), The PyHST2 hybrid distributed code for high speed tomographic
reconstruction with iterative re-construction and a priori knowledge
capabilities. Nucl. Instrum. Methods Phys. Res., Sect. B, Beam Interact.
Mater. Atoms 324, 41--48.

Mitchell, T.M. and Faulkner, D.R. (2012), Towards quantifying the matrix
permeability of fault damage zones in low porosity rocks, Earth Planet.
Sci. Lett. 339-340, 24-31.

Molli, G. and Heilbronner, R. (1999), Microstructures associated with
static and dynamic recrystallization of Carrara marble (Alpi Apuane, NW
Tuscany, Italy), Geologie en Mijnbouw 78, 119-126.

Moura, A., Lei, X.-L., Nishizawa, O., 2005. Prediction scheme for the
catastrophic failure of highly loaded brittle materials or rocks. J.
Mech. Phys. Solids 53, 2435--2455

Nemat-Nasser, S. and Horii, H. (1982), Compression-induced nonplanar
crack extension with application to splitting, exfoliation, and
rockburst, J. Geophys. Res. 87, 6805-6821.

Nemat-Nasser, S. and Obata, M. (1988), A microcrack model of dilatancy
in brittle materials, J. Appl. Mech. 55, 24-35.

Nicksiar, M. and Martin, C.D. (2013), Crack initiation stress in low
porosity crystalline and sedimentary rocks, Eng. Geol. 154, 64-76.

Nicksiar, M. and Martin, C.D. (2014), Factors affecting crack initiation
in low porosity crystalline rocks, Rock Mech. Rock Eng. 47, 1165-1181.

Nozieres, P. (1992), Shape and growth of crystals, in: Godreche, C.
(ed.), Solids Far From Equilibrium, Cambridge University Press,
Cambridge UK.

Odling, N.W.A., Elphick, S.C., Meredith, P., Main, I. and Ngwenya, B.T.
(2007), Laboratory measurement of hydrodynamic saline dispersion within
a micro-fracture network induced in granite, Earth Planet. Sci. Lett.
260, 407-418.

Oesterling, N. (2004), Dynamic recrystallization and deformation
mechanisms of naturally deformed Carrara marble: a study on one- and
two-phase carbonate rocks, PhD thesis, University of Basel.

Ojala I.O., Main, I.G. and Ngwenya, B.T. (2004), Strain rate and
temperature dependence of Omori law scaling constants of AE data:
implications for earthquake foreshock-aftershock sequences, Geophys.
Res. Lett. 31, L24617.

Ollion, J., Cochennec, J., Loll, F., Escude C. and Boudier, T. (2013),
TANGO: A Generic Tool for High-throughput 3D Image Analysis for Studying
Nuclear Organization, Bioinformatics 29, 1840-1841.

Olson, J.E., \& Pollard, D.D. (1991), The initiation and growth of en
echelon veins,~J. Struct. Geol. 13, 595-608.

Ouillon, G. and Sornette, D. (2000), The concept of `critical
earthquakes' applied to main rockbursts with time-to-failure analysis,
Geophys. J. Int. 143, 454--468.

Pal, G., Janosi, Z., Kun, F. and Main, I.G. (2016). Fragmentation by
slow compression of porous rocks, Physical Review E 94, 053003.

Pisarenko, V.F. (1991), Statistical evaluation of maximum possible
earthquakes, Phys. Solid Earth 27, 757-763.

Potyondy, D.O. and Cundall, P.A. (2004), A bonded particle model for
rock, Int. J. Rock Mech. Min. Sci. 41, 1329-1364.

Ramez, M.R.H. and Murrell, S.A.F. (1964), A petro-fabric analysis of
Carrara marble. Int. J. Rock Mech. Min. Sci. Geomech. Abstr. 1, 217-229.

Reches, Z. and Lockner, D.A. (1994), Nucleation and growth of faults in
brittle rocks, J. Geophys. Res. 99, 18159-18173.

Renard, F., Cordonnier, B., Dysthe, D.K., Boller, E., Tafforeau, P. and
Rack A. (2016), A deformation rig for synchrotron microtomography
studies of geomaterials under conditions down to 10 km depth in the
Earth, J. Synchrotron Rad. 23, 1030-1034.

Renard, F., Cordonnier, B., Kobchenko, M. and Kandula, N. (2017),
Microscale characterization of rupture nucleation unravels precursors to
faulting in rocks, Earth Plan. Sci. Lett. 476, 69-78.

Renard, F., Weiss, J., Mathiesen, J. Ben-Zion, Y., Kandula, N. and
Cordonnier, B. (2018), Critical evolution of damage toward system-size
failure in crystalline rock, J. Geophys. Res. 123, 2017JB014964.

Richter, D. and Simmons G. (1974), Thermal expansion behaviour of
igneous rocks, Int. J. Rock Mech. Min. Sci. Geomech. Abstr. 15,
145--148.

Roberts, N.S., Bell, A.F. and Main, I.G. (2015), Are volcanic
seismic\(b\)-values high, and if so when? J. Volc.
Geotherm. Res. 308, 127-141.

Robertson, M.C., Sammis, C.G., Sahimi, M. and Martin, A.J. (1995),
Fractal analysis of three-dimensional spatial distributions of
earthquakes with a percolation interpretation, J. Geophys. Res. 100,
609-620.

Rundle, J.B. and Klein, W. (1989), Nonclassical nucleation and growth of
cohesive tensile cracks, Phys. Rev. Lett. 63, 171-174.

Russ, J.C. (1994), Fractal Surfaces, Plenum Press, New York.

Rydelek, P.A. and Sacks, I.S. (1989), Testing the completeness of
earthquake catalogues and hypothesis of self-similarity, Nature 337,
251--253.

Sammis, C.G. and Ashby, M.F., (1986), The failure of brittle porous
solids under compressive stress states, Actaa Metall. 34 511-526.

Sammis, C.G. and Sornette, D. (2002), Positive feedback, memory and the
predictability of earthquakes, PNAS 99, 2501-2508.

Sammonds, P.R., Meredith, P.G. and Main, I.G. (1992), Role of pore
fluids in the generation of seismic precursors to shear fracture, Nature
359, 228-230.

Schindelin, J., Arganda-Carreras, I., Frise, E., Kaynig, V., Longair,
M., Pietzsch, T., Preibisch, S., Rueden, C., Saalfeld, S., Schmid, B.,
Tinevez, J.Y., White, D.J., Hartenstein, V., Eliceiri, K., Tomancak, P.,
Cardona, A., (2012), Fiji: an open-source platform for biological-image
analysis, Nature Methods 9, 676--682.

Schmittbuhl, J. and Maloy, K.J. (1997), Direct observation of a
self-affine crack propagation, Phys. Rev. Lett. 78, 3888-3891.

Schmittbuhl, J., Schmitt, F., and Scholz, C. (1995), Scaling invariance
of crack surfaces, J. Geophys. Res. 100, 5953-5973.

Sethna, J.P. (2006). Statistical mechanics: entropy, order parameters,
and complexity. Oxford University Press, Oxford UK.

Simmons, G. and Cooper, H.W. (1978), Thermal cycling cracking in three
igneous rocks, Int. J. Rock Mech. Min. Sci. Geomech. Abstr. 15, 145-148.

Smith, R., Sammonds, P.R. and Kilburn, C.R.J. (2009), Fracturing of
volcanic systems: experimental insights into pre-eruptive conditions,
Earth Planet. Sci. Lett. 280, 211-219.

Sornette, D. (2006), Critical phenomena in natural sciences: chaos,
fractals, self-organization and disorder: concepts and tools, Springer,
Berlin.

Sornette, D. (2009), Dragon-kings, black swans and the prediction of
crises, in press in the Int. J. of Terraspace Sci. and Eng. 2, 1-18,
https://arxiv.org/abs/0907.4290

Sornette, D. and Sammis, C.G. (1995), Complex critical exponents from
renormalization group theory of earthquakes: implications for earthquake
predictions, J. Phys. I {[}French{]} 5, 607-619.

Stanley, H.G. (1971), Introduction to phase transitions and critical
phenomena, Oxford University Press, Oxford UK.

Sykes, L.R. and Jaume, S. (1990), Seismic activity on neighboring faults
as a long-term precursor to large earthquakes in the San Francisco Bay
area, Nature 348, 595--599.

Tapponnier, P. and Brace, W.F. (1976), Development of stress-induced
microcracks in Westerly granite, Int. J. Rock Mech. Min. Sci. Geomech.
Abstr. 13, 103-112.

Thomas, A.L., \& Pollard, D.D. (1993), The geometry of echelon fractures
in rock: implications from laboratory and numerical experiments, J.
Struct. Geol. 15, 323-334.

Thouless (1989), Condensed matter physics in less than three dimensions,
in: Davies, P. (ed.), The New Physics, Cambridge University Press,
Cambridge UK.

Turcotte, D. (1997), Fractals and chaos in geology and geophysics,
Cambridge University Press, Cambridge UK.

Tyupkin, Y.S. and Giovambattista, R.D. (2005), Correlation length as an
indicator of critical point behavior prior to a large earthquake, Earth
Planet. Sci. Lett. 230, 85--96.

Utsu, T., Ogata, Y. and Matsu'ura, R.S. (1995), The centenary of the
Omori formula for a decay law of aftershock activity, J. Phys. Earth 43,
1-33.

Vasseur, J., Wadsworth, F.B., Heap, M.J., Main, I.G., Lavallee, Y. and
Dingwell, D. B. (2017), Does an inter-flaw length control the accuracy
of rupture forecasting in geological materials? Earth Plan. Sci. Lett.
475, 181-189.

Vasseur, J., Wadsworth, F.B., Lavallee, Y., Bell, A.F., Main, I.G. and
Dingwell, D. B. (2015), Heterogeneity: the key to failure forecasting,
Sci. Rep. 5, 13259.

Vinciguerra, S., Trovato, C., Meredith, P.G. and Benson P.M. (2005),
Relating seismic velocities, thermal cracking and permeability in Mt.
Etna and Iceland basalts, Int. J. Rock Mech. Min. Sci. 42, 900--910.

Voight, B. (1988), A method for prediction of volcanic eruptions, Nature
332, 125-130.

Voorn, M., Exner, U. and Rath, A. (2013), Multiscale Hessian fracture
filtering for the enhancement and segmentation of narrow fractures in 3D
image data, Comp. Geosci. 57, 44-53.

Weiss, J. (2001), Self-affinity of fracture surfaces and implications on
a possible size effect on fracture energy, Int. J. Fract. 109, 365-381.

Woessner, J. and Wiemer, S. (2005), Assessing the quality of earthquake
catalogues: estimating the magnitude of completeness and its
uncertainty, Bull. Seismol. Soc. Am. 95, 684-698.

Zhang, S., Cox, S.F., and Paterson, M.S. (1994), The influence of room
temperature deformation on porosity and permeability in calcite
aggregates. J. Geophys. Res. Solid Earth, 99, 15761-15775.

Zoller, G., Hainzl, S. and Kurths, J. (2001), Observation of growing
correlation length as an indicator for critical point behavior prior to
large earthquakes, J. Geophys. Res. 106, 2167-2175.

\selectlanguage{english}
\FloatBarrier
\end{document}

