# Deep-Learning Analysis of Fracture Networks Leading to System-Size Failure in Rocks

Joachim Mathiesen, Jessica Mcbeck, Benoît Cordonnier, Yehuda Ben-Zion,

François Renard

► To cite this version:

Joachim Mathiesen, Jessica Mcbeck, Benoît Cordonnier, Yehuda Ben-Zion, François Renard. Deep-Learning Analysis of Fracture Networks Leading to System-Size Failure in Rocks. Journal of Geophysical Research. Machine Learning and Computation, 2024, 1 (3), pp.e2023JH000103. 10.1029/2023JH000103 . hal-05099488

HAL Id: hal-05099488

https://hal.science/hal-05099488v1

Submitted on 5 Jun 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
RESEARCH ARTICLE

10.1029/2023JH000103

**Key Points:**

* Convolutional neural network models predict the proximity of rock samples to failure

* The proximity to failure is predicted better with images than with only fracture density

* Deep learning results find optimal fault angles in the range 10°–30° with respect to $σ_1$, in agreement with empirical failure criteria

**Correspondence to:**

J. Mathiesen and F. Renard,
mathies@nbi.ku.dk;
francois.renard@geo.uio.no

**Citation:**

Mathiesen, J., McBeck, J., Cordonnier, B., Ben-Zion, Y., & Renard, F. (2024). Deep-learning analysis of fracture networks leading to system-size failure in rocks. *Journal of Geophysical Research: Machine Learning and Computation, 1*, e2023JH000103. https://doi.org/10.1029/2023JH000103

Received 12 DEC 2023
Accepted 2 JUL 2024

**Author Contributions:**

**Conceptualization:** Joachim Mathiesen,
Jessica McBeck, Yehuda Ben-Zion,
François Renard

**Data curation:** Benoît Cordonnier

**Formal analysis:** Joachim Mathiesen,
Jessica McBeck

**Funding acquisition:** Joachim Mathiesen,
François Renard

**Investigation:** Joachim Mathiesen,
Jessica McBeck, Benoît Cordonnier,
François Renard

**Methodology:** Joachim Mathiesen,
Jessica McBeck, Yehuda Ben-Zion,
François Renard

**Visualization:** Joachim Mathiesen

**Writing – original draft:**
Joachim Mathiesen, Jessica McBeck,
Yehuda Ben-Zion, François Renard

© 2024 The Author(s). Journal of Geophysical Research: Machine Learning and Computation published by Wiley Periodicals LLC on behalf of American Geophysical Union. This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

<figure><img src="image_2.png" /></figure>

# Deep-Learning Analysis of Fracture Networks Leading to System-Size Failure in Rocks

Joachim Mathiesen¹,² ID, Jessica McBeck² ID, Benoît Cordonnier³, Yehuda Ben-Zion⁴ ID, and François Renard²,⁵ ID

¹Niels Bohr Institute, University of Copenhagen, Copenhagen, Denmark, ²Departments of Geosciences and Physics, The Njord Centre, University of Oslo, Oslo, Norway, ³ESRF, The European Synchrotron, Grenoble, France, ⁴Department of Earth Sciences and Statewide California Earthquake Center, University of Southern California, Los Angeles, CA, USA, ⁵ISTerre, CNRS, IRD, University Grenoble Alpes, University Savoie Mont Blanc, University Gustave Eiffel, Grenoble, France

**Abstract** Fracture networks in rocks and other geomaterials form by seismic and aseismic damage processes that could lead to system-size failure. Here, we use a multi-view convolutional neural network model to predict the stress proximity to macroscopic failure of experimentally deformed rock samples using two-dimensional images. Models are trained on time series of fractures observed in rock samples through dynamic in situ synchrotron X-ray tomography experiments. The results demonstrate that deep learning models outperform traditional estimates based on fracture density, increasing the accuracy of the predictions. Furthermore, the models provide insights into fundamental characteristics of fracture patterns that may provide precursory information on impending material failure. The trained deep learning models estimate the angle of the fracture plane relative to the principal loading direction, which is a key factor contributing to shear failure. The predicted angle of the fracture plane, in the range of 10°–30° with respect to the direction of maximum compressive stress, is consistent with the established failure criteria used in rock mechanics.

**Plain Language Summary** When a rock is subjected to large stress, it can develop fractures that grow and connect to form a network. This network of fractures spans the entire rock right before it fails. Despite extensive research on rock damage, we are still far from being able to predict when a rock will fail. In our manuscript, we use a multi-view convolutional neural network model to identify characteristics of a failing rock. We train our neural network model on images of rock samples exposed to different levels of stress and use it to predict how close a rock is to failure. Our neural network model outperforms traditional estimates based on fracture density. More importantly, our model provides fundamental insight that may offer precursory information on the imminent material failure. For instance, according to our model, the angle of the fracture plane relative to the principal loading direction becomes a key factor contributing to failure.

## 1. Introduction

Laboratory observations of the bulk shape change of cylinders of rock during triaxial compression experiments under conditions of deformation in the brittle crust have highlighted the importance of fracture development on the route to macroscopic system-size failure (Paterson & Wong, 2005; C. H. Scholz, 2019). Since the 1960s and the development of acoustic emission sensors, a large number of experimental studies have shown an increase in acoustic emission activity when approaching system size failure (C. Scholz, 1968; Paterson & Wong, 2005). More recently, triaxial compression experiments with X-ray acquisition at stress and temperature conditions of the upper crust have shown that the total volume of fractures in rocks increases with accelerating rates during increasing differential stress toward macroscopic failure, when the rock experiences a sudden stress drop (Cartwright-Taylor et al., 2020; Renard et al., 2018). When both acoustic and X-ray imaging techniques are combined, the approach to failure shows an increase of seismic emissions and fracture damage and allows quantifying the partition between seismic and aseismic damage (Cartwright-Taylor et al., 2022).

Similarly, several studies have observed variations in geophysical properties and seismicity that precede some large earthquakes that are indicative of an increasing volume of fractures, including changes in water levels in wells, variations in seismic wave speeds and attenuation, and localization of low-magnitude earthquakes and rock damage (Ben-Zion & Zaliapin, 2019, 2020; Bufe & Varnes, 1993; Roeloffs, 1988; C. H. Scholz et al., 1973).
**Writing – review & editing:**
Joachim Mathiesen, Jessica McBeck,
Yehuda Ben-Zion, François Renard

These changes may arise before large earthquakes because small earthquakes or aseismic slips or both produce
increasing crack density and permeability surrounding the large rupture zone as well as within it.

The clear relationship between fracture development and macroscopic failure in laboratory experiments (e.g., Paterson & Wong, 2005; C. H. Scholz et al., 1973), and more tenuous relationship in crustal faults (Ben-Zion, 2008; C. H. Scholz et al., 1973), have led some studies to suggest that a critical volume of fractures is necessary to trigger macroscopic failure (Brantut et al., 2013; Kachanov & Sevostianov, 2012; Kranz & Scholz, 1977; Lockner et al., 1991). Furthermore, the spatial organization of these fractures, such as their clustering or localization and their orientation relative to the principal stresses, likely influences the ability of a particular fracture network to trigger macroscopic failure (Ben-Zion & Zaliapin, 2020; Kachanov & Sevostianov, 2012).

Following these ideas, previous work has used time series of three-dimensional fracture networks observed in X-ray tomography experiments and machine learning to predict the proximity to macroscopic system-size failure in rock core samples (J. A. McBeck et al., 2020). This approach requires calculating specific geometric properties for all of the fractures at all of the differential stress steps (and tomograms) of each experiment. The choice of geometrical properties to analyze is based on pre-existing failure criteria and intuitions about which characteristics may be the most helpful for predicting the proximity to failure. Specific geometric properties may include, for example, fracture volume, dimensions and orientation, distance to the nearest fracture, and the Gini coefficient of the fracture population. If some important and unknown properties that control system size failure are not selected, this type of supervised analysis can unintentionally ignore important characteristics of the fracture network that signal the proximity to macroscopic failure.

Several recent studies used machine learning techniques to analyze data recorded in laboratory experiments. Jasperson et al. (2021) and Shreedharan et al. (2021) demonstrated success in predicting large stick-slip events in laboratory frictional experiments with deep learning of acoustic emission data and estimated changes of elastic properties recorded in the experiments. Wang et al. (2022) and Laurenti et al. (2022) used transformer networks to forecast evolving stress and displacement in laboratory frictional experiments. Karimpouli et al. (2024) used unsupervised learning and features extracted from acoustic emission catalogs to forecast large slip events in triaxial stick-slip experiments. In this study, we continue these efforts to predict the stress proximity to macroscopic failure in laboratory fracturing experiments by using a convolutional neural network to analyze two-dimensional slices of three-dimensional segmented tomograms of eight experiments. This analysis does not require identifying the characteristics of the fracture network that may be useful for predicting the proximity to failure a priori, nor does it require simplifying the complex spatial distribution of fractures that differ in volume by orders of magnitude (Cartwright-Taylor et al., 2020; Kandula et al., 2019; Renard et al., 2018) to a few basic characteristics.

The developed deep learning model is shown to provide predictive information on the proximity to failure based on sets of two-dimensional images of fracture patterns where the fractures have been extracted by a pre-processing segmentation step. Moreover, the model predicts that rock samples that contain fractures oriented at 10°–30° with respect to the direction of the main compressive stress have a higher proximity to failure, in agreement with existing failure criteria used for rocks under compression (e.g., Anderson theory of faulting Anderson (1905), see also Paterson and Wong (2005)).

## 2. Data, Method, and Model Training

Our initial data were taken from a series of rock deformation experiments where synchrotron micro-computed tomography data were acquired during the approach of rock samples to failure. Eight rock core samples were deformed within the Hades triaxial rock deformation apparatus until system-size failure occurred. This apparatus, installed on the beamline ID19 at the European Synchrotron Radiation Facility, enables the performance of rock deformation experiments under in situ conditions of pressure and temperature, as well as the imaging of rock samples using time-lapse X-ray microtomography (Renard et al., 2016). The eight samples were described in previous studies and the data are publicly available (Kandula et al., 2019; J. A. McBeck et al., 2020; Renard et al., 2019, 2018). The eight core samples are two Carrara marbles (marble 1, 2), three monzonites (monzonite 3, 4, 5), and three Westerly granites (granite 1, 2, 4), see (J. A. McBeck et al., 2020). All are crystalline rocks with porosity on the order of or below 1%. The core samples were 10 mm height, with a diameter of either 5 mm (marble) or 4 mm (granite, monzonite).
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Sketch visualizing how we extract either one or two perpendicular two-dimensional slices from a three-dimensional binarized tomogram where fractures have been extracted by a segmentation procedure. We then develop two convolutional neural networks trained with either one or two slices, and we call the latter a multi-view approach. The core sample is a rock cylinder and the fractures are represented in gray (left image) or black (middle and right image).</figcaption></figure>

First, we summarize the experimental protocol. Experiments were carried out at room temperature (22–24°C) and each rock core sample was surrounded by a jacket on which a constant confining pressure was applied in the range 5–35 MPa with oil. No pore fluid pressure was applied. The axial stress, $σ_1$, on each sample was increased by steps until system-size failure occurred. After each stress step increase, a three-dimensional X-ray acquisition was performed under constant stress, and the entire rock sample was imaged to produce a tomogram with a spatial sampling distance of 6.5 µm, corresponding to the voxel size, and a spatial resolution equal to two to three times this sampling distance. The last tomogram was acquired within 1 MPa before the axial stress at failure. Tomograms have dimensions in the range 2,000 × 1000 × 1,000 voxels. After tomography reconstruction of the samples, each voxel records X-ray adsorption coded with a 16-bit grayscale value, with dark colors corresponding to voids (e.g., fractures, pores) and lighter colors corresponding to the minerals of the rocks. We use the procedure described in (J. A. McBeck et al., 2021) to extract microfractures (i.e., segmentation) and transform the tomograms into binary volumes where the value of one corresponds to a void (e.g., fracture or initial pore) and zero to a mineral (see Figure 1). The stress-strain curves of the eight experiments, the confining pressure, and the number of tomograms used are shown in Figure 2.

In principle, it would be possible to train a model on the full binary three-dimensional data. However, the performance of the model would suffer from lack of data due to the sparseness of fractures, which amount to less than a few percent of the rock volume, even when approaching failure. Various suggestions have been proposed on how to draw inferences based on three-dimensional volumetric shapes (Qi et al., 2016; Su et al., 2015). Inspired by the work of Su et al. (2015) who proposed a technique to represent three-dimensional data by using two-dimensional images, and that these authors called a multiple views approach, we consider two models (see Figure 1). The first is trained on individual two-dimensional slices of the three-dimensional tomograms, whereas the second model adopts the multi-view approach and is trained simultaneously on two two-dimensional slices cut at different angles. For simplicity, we limit ourselves to slices that are perpendicular and aligned with the principal direction of compression, $σ_1$. Note that in contrast to most multi-view approaches, we examine data on slices cut through the sample instead of observing the fracture network from the outside.

The two-dimensional images used in the training have been produced in the following way. For each sample, we compute the center of mass of the fracture volume in the last stress step before failure. We designate the direction along the principal loading, $σ_1$, as the z-axis. We cut two-dimensional slices with normal vectors in the x-y plane using the center of mass as the center of coordinates. The two-dimensional slices are then down-sampled to approximately half size using a density-preserving Gaussian filter, resulting in grayscale images of 600 by 250 pixels. We could train the model on the full-scale images, but again, due to the sparseness of the fracture volume little information is lost by down-sampling and at the same time the training performance is greatly improved. Thus, from an initial image that contains only integer numbers 0 and 1, the pixels now have real numbers. In
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Stress-strain curves of the eight rock samples in the final stages before failure, within the last 20% of the stress at failure and that we use in the present study. Samples were deformed in the Hades rig at the European Synchrotron Radiation Facility. The confining pressure of each experiment, $\sigma_2$, is indicated. The black circles indicate the acquisition of a three-dimensional tomogram. Labels of the x- and y-axis indicated on the top left plot apply to all the other plots. The complete stress-strain curves are shown in the Figure A1 in the Appendix A.</figcaption></figure>

Appendix C, we show that the choice of vertical slices has only a minor effect on the capacity of the model to detect shear planes oriented at an angle $\theta \sim 30^{\circ}$ to the direction of compression.

The target value that we would like to predict with our model is the proximity to the stress at failure. The same approach was used by C. Scholz (1968), who defined the proximity to failure by the percentage of fracture stress. In the following, we define a slightly different empirical non-dimensional number:

$$ \Delta_t = \frac{\sigma_f - \sigma}{\sigma_f}. \qquad (1) $$

Here we use the value of the differential axial stress applied on each sample at a given stress step, $\sigma$, and the stress applied at the moment when a sample in the experiment fails, $\sigma_f$. Thus, $\Delta_t = 0$ corresponds to failure and $\Delta_t = 1$ indicates the start of the experiment where the axial differential stress is zero. We have considered values of $\Delta_t$ in the interval [0 0.7]. However, because the fracture density is extremely low for $\Delta_t$ larger than 0.2, we focus here on the interval $\Delta_t = [0 0.2]$, corresponding to the last 20% of the stress-strain curve shown in Figure 2. In this interval of target values, we end up with a data set consisting of a total of 40,000 images distributed across the three rock types granite, marble, and monzonite. The full stress-strain curves are included in the Figure A1 in Appendix A.

The models we use are, respectively, a standard convolutional neural network applied to individual slices and a multi-view convolutional neural network applied to two perpendicular slices simultaneously. In Figure 3, we illustrate the multi-view model. The standard convolutional neural network has a similar structure except that the module accounting for the second input is deactivated during the training. Both convolutional neural networks contain five convolution layers and end with a three-layered densely connected output module with a final output layer using a rectified linear activation unit (ReLU); see (Figure 3). The ReLU activation ensures that our model does not predict negative values for $\Delta_p$. The individual convolution layers all have ReLU activations and consist of 64 filters of size 3 × 3 with a unit stride. We use padding to maintain the dimension of our images after each convolution. The convolutional layers are each followed by a max pooling layer. The pooling layers similarly use a unit stride. The models have been trained with batch-normalization before the non-linearity in each layer.

$$ \text{RMSE} = \sigma_{sd}(\Delta_p - \Delta_t), \qquad (2) $$
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Network topology used in the deep learning model to predict proximity to failure from respectively one and two slices. When using one slice, only one of the two convolutional neural network models is used. The network flattens the convolutional layers, which feed into a dense neural network used for the final prediction.</figcaption></figure>

where $\Delta_p$ is the predicted proximity to failure and the function $\sigma_{sd}$ is the standard deviation taken over all images in either the training or the validation pool, respectively.

In the following, we present results by first training our model across all rock materials, Approach 1 (abbreviated as **Ap1** in the following), and afterward by training on data where we exclude in turn the different samples, Approach 2 (**Ap2**). The latter approach offers a more conservative estimate of the performance to be expected from neural networks and in addition provides realistic information about the generality of predicting fractures across materials types. Results from **Ap1**, the distribution of $\Delta_p$, is presented in Figures 4 and 5 summarizes all our results and compares **Ap1** and **Ap2**.

In **Ap1**, we pool all data across the rock samples and divide the data in two groups the training data (70%) and the final test data (30%). We reserve 15% of the training data for validation (which therefore does not contribute to the training of model parameters). In approach **Ap2**, we train eight models by in turn excluding one of the eight data sets during training. The excluded data is only used in the final testing of our model. Among the slices in the seven samples used for training, we divide again the 85% in a training set and the remaining 15% in a validation set.

Although trivial, we emphasize that before each training run, our model is reset, so that the model contains no information about the rock sample used for the final testing.

We construct and train our models using TensorFlow through the KERAS API (Chollet, 2015). For the training, we have used both the Nesterov-accelerated Adaptive Moment (NADAM) and the root mean square propagation (RMSProp) optimizers. The RMSProp provided overall more consistent results whereas NADAM often did not converge very well to a good model. In the Appendix B (Figure B1), we show the optimal training history used for the model in **Ap1**. Our models are generally trained by up to approximately one hundred epochs using adaptive learning rates that decrease as the network training converges. Training was stopped when there was negligible improvement in the overall accuracy of the model (less than 5% on average over 10 epochs). In general our training runs all stopped between 40 and 100 epochs for **Ap1**, where the training stopped 10–30 epochs for **Ap2**. In general, we did not see any improvement for **Ap1** by going beyond 100 epochs. Whereas the training loss continues to improve, the validation loss quickly reaches a plateau and continues to fluctuate and only slowly
<figure><img src="image_6.png" /><figcaption>Figure 4. Histograms of the difference between the predicted proximity to failure, Δp, and the true value of this parameter, Δt in Ap1. We compare the histogram of the two-slice model (CNN Double) and the prediction based on the fracture density alone. The results are achieved from models trained across all material types and samples.</figcaption></figure>

starts to increase in value. We believe that because of the simplicity and to some extent the generality of the fracture patterns, the models are not quickly ending in a state with significant overfitting. The learning rate used in the models was initially set to 0.001. We introduce an automatic reduction by a factor 0.9 in the learning rate when a plateau in the training loss is reached (i.e. no improvement over five epochs, the patience parameter set to 5). The reduction in learning rate only had a minor influence on the convergence during training of Ap1 and did not influence the training of Ap2. The training of the models was performed on a multi-GPU computer using a batch size of 256 images for both the training and validation data.

### 3. Results

The main result is presented in Figure 5, where we compare the predicted proximity to failure with the true value Δt for both Ap1 and Ap2.

As a first test of the predictive power of our models, we compare the performance with a baseline model based on random guesses of the proximity to failure. The random guessing is made in the following way. We take all experimental loading steps that are within 20% of failure, which means $0.0 \le \Delta_t \le 0.2$. For each slice from a given experiment, we pick, among all the loading steps, a random value of the proximity to failure Δp. The true proximity to failure is given by Δt. The root mean square error (RMSE) of our random guess is then found by averaging over all slices in an experiment

$$ \sqrt{\langle (\Delta_t - \Delta_p)^2 \rangle} = \sqrt{\langle \Delta_t^2 \rangle + \langle \Delta_p^2 \rangle - 2\langle \Delta_t \rangle \langle \Delta_p \rangle} = \sqrt{2\langle \Delta_t^2 \rangle - 2\langle \Delta_t \rangle^2} \quad (3) $$

Here, the random guesses are taken from the population of true values Δt; therefore, the means of the true and guessed values are equal, $\langle \Delta_t \rangle = \langle \Delta_p \rangle$, as well as the second moments. If the loading steps in the experiment were

uniformly distributed, we would find that $\sqrt{\langle (\Delta_t - \Delta_p)^2 \rangle} = 2\sqrt{1/150 - 1/200} \approx 0.082$. As observed in Figure 5, the null-model (black squares) is approximately close 0.082 for all the experiments (the loading steps are fairly uniformly distributed).

As a second test of the predictive power of neural network models, we use the fracture density within each experiment to estimate the proximity to failure. For each experiment, we compute the median fracture density for each loading step (within 20% of failure) and produce the corresponding curve plotting the step-wise fracture density versus the proximity to failure. For each slice in the same experiment, we then find the point on this curve that has the closest fracture density and read off the proximity to failure, which is the estimated proximity to failure. We plot this estimate in Figure 5. We finally compare our single slice convolutional neural network and the multi-view convolutional neural network (for Ap1 and Ap2) with our null-model and the density estimate.

In Ap2, we find that for all but one experiment (granite 2), the convolutional neural networks perform better than baseline models. The multi-view convolutional neural network performs better for all experiments than the baseline models and Ap1 greatly outperforms Ap2. The latter is to be expected since we strongly restrict knowledge about the specific material during training in Ap2. The reason that Ap2 performs worse than the null-model for granite 2, is likely due sub-optimal training (staying faithful to our stopping criterion) and due to the relatively few loading steps for this specific material.

<figure><img src="image_7.png" /><figcaption>Figure 5. Root-mean-square of the difference between the predicted proximity to failure, Δp, and the true value of this parameter. We present results from four classes of models: a null model where the proximity to failure is randomly assigned from a uniform distribution, a model based on the density of fractures when approaching failure, and the two convolutional neural networks with one (CNN Single) or two (CNN Double) slices. The figure contains results from CNN Double in both our approaches Ap1 and Ap2.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Predicted proximity to failure using synthetic images. We generate images of randomly positioned micro-fractures of one pixel size. The resulting speckle noise image looks like the 600 × 250 pixels image shown in the right panel. In the left panel, we show in different colors the true value $\Delta_t$, as function of the fracture density for the eight samples used to train our model. The black line is the predicted proximity to failure achieved by simply fitting our **Ap1** model to the synthetic images. In general, the prediction based on the synthetic images follows the trend of the experimental data, getting closer to failure with increasing fracture density. In other words, density, not unexpectedly, does contain significant information on the proximity to failure.</figcaption></figure>

In the Figures 6 and 7, we explore the capabilities of our fully trained network model by applying it to synthetic images generated in two different ways. We use the synthetic images to identify characteristic fracture patterns that our will be more sensitive to than others. In Figure 6, we generate images consisting of small fractures introduced throughout the image. These synthetic fractures are selected at random locations and have a length of one pixel. For different fracture densities or non-zero pixel densities, we use our model trained in **Ap1** to predict the supposed proximity to failure for the synthetic images. The predicted values $\Delta_p$ are represented by the circles connected by the black line. The value of each black circle in Figure 6 is derived from an average of predictions over one hundred synthetic images. Overall, we observe that an increase in fracture density results in a prediction of a closer proximity to failure but we also observe that this effect is not very strong for intermediate fracture densities. The colored points in the background of the figure represent the proximity to failure versus fracture density for the real rock data. The colors represent the individual rock samples (see the figure legend).

In Figure 7 we consider another set of synthetic fracture images. Whereas we considered micro-cracks of unit length in the previous example, we shall here consider fractures extending to the system size in the form of stripe patterns that are orientated at an angle $\theta$ with the principal stress direction. The results on the synthetic images (Figure 7) show that the stripe patterns (i.e., fractures) oriented at 10°–30° to $\sigma_1$ are predicted to be at closer proximity to failure. Furthermore, the density of the stripes (i.e., the fracture density) has relatively little effect, similar to the observation for the intermediate fracture densities in Figure 6.

To verify that fracture orientation does have an effect on the distance to failure, we have extracted all fractures from the 2D slices and calculated the orientation of their main elongation axis using a principal component analysis (Figure 9). Results show that, when approaching failure, the distribution of orientations accumulates significantly in the range 0°–20°, in agreement with the results shown in Figure 6.

Finally, in Figure 8 we present a few examples of patterns that trigger the filters in our model. The patterns are generated in the following way. We start from images of random white noise and apply a gradient ascent function,
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Predicted proximity to failure, $\Delta_p$, as function of fracture orientation. Based on a model trained on single slices from all rock types, results show the fracture plane orientation relative to the direction of principal stress that is closest to failure. The predictions are computed by presenting to our model fracture patterns (600 pixels by 250 pixels, shown in the right panel) with different orientations, $\theta$. The different colors in the left panel correspond to fracture patterns with different spacing ($\Delta w$) between the stripes. Results show that optimal fracture orientation is within the range 10°–30° from the main compressive stress, in agreement with a wide range of experimental rock physics studies (e.g., McCormick & Rutter, 2022) and the empirical Coulomb-Anderson failure criterion (Anderson, 1905).</figcaption></figure>

which iteratively shifts the white noise toward a pattern that maximizes the activation of individual filters. In the figure, we show six randomly selected filters from the top layer in our model. One of the filters seem to pick up a strong sense of directionality (the upper right panel in Figure 8). The directional pattern is consistent with our prediction of the optimal fracture angle when close to macroscopic failure; see Figure 7.

# 4. Discussion

## 4.1. Machine Learning and Laboratory Earthquakes

In recent years, several studies have developed machine learning models to predict the occurrence of small failure events in laboratory experiments. Using acoustic emissions recorded in friction experiments where a layer of gouge is sheared between two blocks, a study applied a random forest model to predict the time to the next failure event (Rouet-Leduc et al., 2017). The method was extended to detect slow and fast failures using combinations of supervised and unsupervised learning methods (Hulbert et al., 2019; Johnson et al., 2021). A deep learning technique that uses acoustic emission data as input to predict the shear stress on the fault and the time to failure outperforms other machine learning techniques, such as decision tree methods or gradient boosting algorithms (Laurenti et al., 2022). In an analog experiment of a subduction zone (Corbi et al., 2019), a gradient boosting method was developed to predict the proximity to the next failure events using surface displacement data.

In the present study, we apply convolutional neural network models to predict the proximity to failure using two-dimensional images of fracture networks in experiments where rock samples are loaded until system-size failure. Our approach has the advantage that we do not have to make prior assumptions about the relevant features to include in the model, as previously done (J. A. McBeck et al., 2020). The relevant features and the
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Visualization of a few selected filters in the second convolutional network layer. We start from a random image with white noise and iteratively apply gradient ascent steps in order to maximize the activation of a single filter in the target layer. In this way, we enhance features in the random image that strongly activate a given filter. Some of the filters respond to stripe patterns with a specific angular orientation, whereas others do not respond to any discernible features.</figcaption></figure>

important physics emerge through the training. For example, instead of measuring individual fracture angles relative to the principal direction of loading and using them as inputs to a model, the importance of the angular fracture orientation when approaching the macroscopic failure emerges during the training of the model (as seen in Figure 7). In other words, and in a broader perspective, the models we use have the potential to uncover physical mechanisms not necessarily anticipated in more traditional approaches, where preselected characteristic variables are identified and used in a model.

## 4.2. The Emergence of Optimal Shear Fracture Angles in Deep Learning Models

Rock physics experiments where samples are loaded under compression show that following the linear stage, at about half of the stress of macroscopic failure, the rock begins to dilate, and this dilation continues until macroscopic failure (Paterson & Wong, 2005). Previous work attributed this dilation to fracture opening that primarily occurs perpendicular to the direction of the main compressive stress, $\sigma_1$, as the fractures lengthen parallel to the $\sigma_1$ direction (Brace et al., 1966). Initially, these fractures have various orientations and then coalesce into a connected fracture network that leads to system-size failure of the rock sample. The angle between this network and $\sigma_1$ is around 20° (Reches

& Lockner, 1994), and can be predicted using bifurcation theory (Rudnicki & Rice, 1975). Recent experiments on samples of Carrara marble and Solnhofen limestone have shown that this angle varies in the range 10°–30° under compression and slightly increases with confining pressure (McCormick & Rutter, 2022). We cannot test the effect of confining pressure using our data. However, a range of optimal angles produces similar values of $\Delta_f$ in synthetic data; in particular, a range of 10°–30° emerges from the models (Figure 7), similar to the range in experiments of McCormick and Rutter (2022).

In a previous study where we used supervised machine learning, we extracted the geometric properties of the fracture network as features to predict the proximity to failure in terms of the differential stress (J. A. McBeck et al., 2020). The results showed that the fracture characteristics that most strongly control the model predictions are the orientation of the fractures relative to the $\sigma_1$ direction, the shape anisotropy or thinness of the fractures, the fracture aperture and the spacing between the fractures. These characteristics control the strain energy density of the system, and so this metric may provide more robust predictions of the proximity to failure than other failure criteria, such as the Mohr-Coulomb criteria. This analysis did not find that the total fracture volume in a tomogram or subvolume helps the models predict the proximity to failure in the rocks with more heterogeneous distributions of fractures, the Westerly granite and quartz monzonite. This analysis included features derived from the distance between the geometric centroids of the fractures, but this metric is limited because it does not consider any information about the size of the fractures. Consequently, the spatial organization of many smaller fractures around one much larger fracture may appear similar in the derived localization metrics as the organization of many smaller fractures around another smaller fracture. Quantifying the location using geometric centroids may also lead to incomplete approximations of the spatial distribution because the tip of a large fracture could be very close to the tip of another fracture, but the centroids of the two fractures could be far apart and thus appear distant to the models. On the contrary, the model we develop here using deep learning does not rely on the choice of features and predicts that the optimal orientation of fracture that controls failure is in the range 10°–30°, in agreement with experimental results (McCormick & Rutter, 2022).

<figure><img src="image_11.png" /><figcaption>**Figure 9.** Probability density function of fracture orientations in the range 0°–90°, where 0 indicates that the fracture is parallel to the loading direction, for three different proximities to failure. When approaching failure, fractures accumulate in the orientation range 0°–20°, in agreement with results shown in Figure 6.</figcaption></figure>
A recent study suggests that the dominant deformation mode immediately preceding macroscopic failure in rock cores of Westerly granite is tensile as a larger volume contains fractures with orientations parallel and subparallel to $σ_1$ than within the range expected for shear (J. McBeck et al., 2023). In our deep-learning model, some synthetic data with a higher stripe density (e.g., higher fracture density) and orientations of 5° produce a closer proximity to failure than synthetic data with other orientations (Figure 7). This result agrees with older experimental studies that have shown that system-size failure occurs by the linkage of tensile fractures (oriented parallel and sub-parallel to $σ_1$) along the future fault plane that develops at 10°–30° with respect to $σ_1$ (Peng & Johnson, 1972).

### 4.3. Using Neural Networks to Estimate Proximity to Failure

In order to evaluate the quality of our neural network models, we compare them up against a null model. The null model that we considered here, is of the most simply type where we randomly guess the proximity to failure. When limiting our training data to samples within the 20% stress at failure, this null model leads to a root-mean-square error of 0.08. When we use fracture density as a predictor, we do not improve much relative to the null-model. We note that the fracture density varies greatly between different rock types and could be to some extent sensitive to the way the 3D-image reconstruction is performed from the raw X-ray data (J. A. McBeck et al., 2021). We considered two approaches for training a convolutional neural network model, one where we do not train on the same material type that we want to draw inference about and another where we include the same material type in the training. The former approach **Ap2** provides a more conservative estimate of the proximity to failure (with an RMSE in the range 0.04–0.09). The latter approach **Ap1** offers an RMSE in the range 0.01–0.03. With more data from additional experiments across many different material types, we could greatly improve performance and probably see the conservative estimate approaching that of **Ap1**.

### 4.4. Future Studies and Relevance to Crustal Data

An obvious and labor-intensive improvement of our model would be to acquire and use more experimental rock samples to train the deep learning model. With more experiments on a wider range of rocks, it could be tested if the deep learning model could predict the rock type based on fracture patterns. Another improvement would be to train the model on smaller three-dimensional sub-volumes instead of two-dimensional slices. This approach may yield better predictions of the correlation between the orientation of the fault plane and the proximity to failure. Finally, with more samples available, we could include rock type as a categorical input and perhaps other variables, such as overall fracture density and confining pressure, to better predict the proximity to failure from two-dimensional images of fracture patterns. It might also be possible to test our model on field-scale data, for example, on fracture networks extracted from optical images of outcrops using machine learning methods (Mattéo et al., 2021).

Extrapolating the results of the present study to larger scales and natural fracture data sets would require more experimental observations to cover a wider range of rock types and stress states, and testing the transfer learning capacity of our model to field data. A consequence of fault networks is that natural seismicity is strongly clustered in the crust (Zaliapin & Ben-Zion, 2022), and often aligns along several dominant directions exhibiting an anisotropic structure (Ross et al., 2022). Moreover, seismicity may occur on faults at depth that do not outcrop at the surface and are thus partially hidden. A recent study has proposed to automatically map fault networks from field images and topographic data using deep learning (Mattéo et al., 2021). The results are maps of outcropping fracture networks that contain a similar information as the images used in the present study: a digitized network of fractures. The comparison of these data sets with a large series of X-ray tomography experiments could be the topic of a future study that measures the proximity to failure of natural fracture sets and the orientation of the main compressive stress that would minimize this proximity. Applying our method on images of faults and fractures at the surface of the earth, in conjunction with information from seismic catalogs and other data sources, could possibly be used to detect partially-hidden large faults. This may be the subject of a follow up work.

## 5. Conclusions

Using fracture patterns extracted from X-ray tomograms of rock core samples approaching failure under confinement, we developed a deep learning approach to estimate the proximity to failure based on two-dimensional images (Figure 1). An important finding is that using images from one or two cross sections is
sufficient to provide an estimate of the proximity to failure. Relying solely on the density of fractures in the models provides less reliable estimates of this proximity than using the fracture patterns present in the two-dimensional images. Additionally, the deep learning model trained on experimental data also reveals that the optimal orientation of fractures leading to failure is in the range of 10°–30°, and in data with the highest fracture density near 5°, in agreement with empirical theories of shear faulting (e.g., Anderson, 1905) and recent laboratory observations (J. McBeck et al., 2023; McCormick & Rutter, 2022). Our deep learning approach reproduces key features of rock mechanics using images as input and without prior knowledge of this theory. Extending this work to larger scale natural data would require performing experiments for a wider range of rocks and stress conditions such that experimental fracture networks may be compared with natural fracture networks that surround fault zones.

## Appendix A: Complete Stress-Strain Curves of the Eight Samples

See Figure A1.

<figure><img src="image_12.png" /><figcaption>**Figure A1.** Complete stress-strain curves of the eight rock samples. The black circles indicate the acquisition of a three-dimensional tomogram. The red circles and inset show the tomograms (number in red) located within 20% of the stress at failure that we used in the present study. These data are also shown in this figure. The labels of the x- and y-axis indicated on the top left plot apply to all the other plots.</figcaption></figure>

## Appendix B: The Fully Trained Model—Losses and Results

The training and validation loss for our model trained across all material types are shown in Figure B1. The detailed results summarized in the histogram in Figure 4 are shown for the eight samples in Figure B2.
<figure><img src="image_13.png" /><figcaption>**Figure B1.** Training and validation loss as function of the number of epochs for the **Ap1** (upper panel) and **Ap2** (lower panel) models, respectively. In the upper panel, the minimum validation loss is reached after 60–70 epochs. After that point the validation is not improved within the training fluctuations, whereas the training loss continues to be reduced. In the lower panel, the models were trained on rock type not used in the validation. We present training data for the three rock types. The training reaches a minimum validation loss after 10–30 epochs. Note that the validation drops from a root mean square error (RMSE) similar to that of the null model and to a RMSE value around 0.04–0.06 as presented in Figure 5.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>**Figure B2.** Individual prediction across the experimental samples used in approach Ap1 in Figure 5.</figcaption></figure>

## Appendix C: Sub-Sampling a Three-Dimensional Volume

In general, sub-sampling a three-dimensional volume by two-dimensional slices would lead to spurious predictions of the fault plane orientation. Following the Mohr-Coulomb criterion, the fault plane developed close to rock failure is expected to form at a critical angle $\theta$ with the largest principal stress direction $\sigma_1$. The planar slices we use in the training are aligned with the principal loading direction. If we define a coordinate system with a z-axis along the principal loading, our planar cuts will have normal vectors confined to the xy-plane,

<figure><img src="image_15.png" /><figcaption>**Figure C1.** Probability density function, panel (a), and cumulative distribution function, panel (b), of the apparent fracture plane angle forming from synthetic random two-dimensional slices along the dominant principal stress direction.</figcaption></figure>
$$ \mathbf{n}_c = \begin{pmatrix} \cos \phi \\ \sin \phi \\ 0 \end{pmatrix} $$

The planar cut with the fracture plane will form a line, which in general will not preserve the angle with the z-axis, that is, the apparent fracture angle $\tilde{\theta}$ will follow from the relation

$$ \cos(\tilde{\theta}) = \mathbf{e}_z \cdot (\mathbf{n}_f \times \mathbf{n}_c) = \cos \theta \sin(\phi) $$

or similarly

$$ \tilde{\theta} = f(\phi, \theta) = \arccos[\cos \theta \sin(\phi)] $$

For a fixed fracture angle $\theta$ we get that random planar cuts (where $\phi$ will be uniformly distributed between 0 and $\pi$) will form apparent fracture angles distributed according to the probability density function

$$ p(\tilde{\theta}) = \frac{1}{\pi} \frac{1}{|\frac{\partial \tilde{\theta}}{\partial \phi}|} = \frac{1}{\pi} \left| \frac{\cos \theta \cos \phi}{\sqrt{1 - \cos^2 \theta \sin^2 \phi}} \right|^{-1} $$

where

$$ \phi = \arcsin\left[\frac{\cos\tilde{\theta}}{\cos\theta}\right] $$

and therefore

$$ p(\tilde{\theta}) = \frac{1}{\pi} \frac{1}{|\frac{\partial \tilde{\theta}}{\partial \phi}|} = \frac{1}{\pi} \sqrt{\frac{1 - \cos^2 \tilde{\theta}}{\cos^2 \theta - \cos^2 \tilde{\theta}}} $$

The distribution is confined to an interval between the faulting angle $\theta$ and $180 - \theta$ and diverges near both end of the interval. The latter angle corresponds to a system turned upside down. We use here that the faulting angle is expected to be close to 30°. The probability density function is shown in Figure C1 together with the cumulative distribution. More than 30% of the data are within 10 degrees of the critical angle. So despite the fact that we throw away important information by training on two-dimensional thin slices, we still recover information about the angle of the micro-fault planes forming since the distribution following from random sampling has a distinct peak around the characteristic angle.

## Conflict of Interest

The authors declare no conflicts of interest relevant to this study.

## Data Availability Statement

The dynamic three-dimensional microtomography data for the eight samples are available as two Carrara marble samples (Renard, 2019), monzonite 3 and monzonite 4 (Renard, 2017), monzonite 5 (Renard, 2018), and three granite samples (Renard, 2021).
Acknowledgments

Beamtime was allocated at the European Synchrotron Radiation Facility (Long Term Proposal ES-1190). This project has received funding from the Research Council of Norway (project "Emergent Networks," Grant 300435), the Norwegian e-infrastructure for Research and Education UNINETT Sigma2 (project NS9073K), and the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement no. 101019628 BREAK). F. Renard and J. McBeck acknowledge support from the project FricFrac funded by the Center for Advanced Study (CAS) at the Norwegian Academy of Science and Letters during the academic year 2023-2024. The work of Y. Ben-Zion on analysis and interpretation was supported by the U. S. Department of Energy (Award DE-SC0016520). J. Mathiesen acknowledges support from the ESS lighthouse on hard materials in 3D, SOLID, funded by the Danish Agency for Science and Higher Education, Grant 8144-00002B.

References

Anderson, E. M. (1905). The dynamics of faulting. *Transactions of the Edinburgh Geological Society*, 8(3), 387–402. https://doi.org/10.1144/ transed.8.3.387

Ben-Zion, Y. (2008). Collective behavior of earthquakes and faults: Continuum-discrete transitions, progressive evolutionary changes, and different dynamic regimes. *Reviews of Geophysics*, 46(4), 4006. https://doi.org/10.1029/2008rg000260

Ben-Zion, Y., & Zaliapin, I. (2019). Spatial variations of rock damage production by earthquakes in southern California. *Earth and Planetary Science Letters*, 512, 184–193. https://doi.org/10.1016/j.epsl.2019.02.006

Ben-Zion, Y., & Zaliapin, I. (2020). Localization and coalescence of seismicity before large earthquakes. *Geophysical Journal International*, 223(1), 561–583. https://doi.org/10.1093/gji/ggaa155

Brace, W., Paulding Jr, B., & Scholz, C. (1966). Dilatancy in the fracture of crystalline rocks. *Journal of Geophysical Research*, 71(16), 3939–3953. https://doi.org/10.1029/jz071i016p03939

Brantut, N., Heap, M., Meredith, P., & Baud, P. (2013). Time-dependent cracking and brittle creep in crustal rocks: A review. *Journal of Structural Geology*, 52, 17–43. https://doi.org/10.1016/j.jsg.2013.03.007

Bufe, C. G., & Varnes, D. J. (1993). Predictive modeling of the seismic cycle of the greater San Francisco bay region. *Journal of Geophysical Research*, 98(B6), 9871–9883. https://doi.org/10.1029/93jb00357

Cartwright-Taylor, A., Main, I. G., Butler, I. B., Fusseis, F., Flynn, M., & King, A. (2020). Catastrophic failure: How and when? Insights from 4-D in situ X-ray microtomography. *Journal of Geophysical Research: Solid Earth*, 125(8), e2020JB019642. https://doi.org/10.1029/ 2020jb019642

Cartwright-Taylor, A., Mangriotis, M.-D., Main, I. G., Butler, I. B., Fusseis, F., Ling, M., et al. (2022). Seismic events miss important kinematically governed grain scale mechanisms during shear failure of porous rock. *Nature Communications*, 13(1), 6169. https://doi.org/10.1038/ s41467-022-33855-z

Chollet, F. (2015). *Keras*. GitHub. Retrieved from https://github.com/fchollet/keras

Corbi, F., Sandri, L., Bedford, J., Funiciello, F., Brizzi, S., Rosenau, M., & Lallemand, S. (2019). Machine learning can predict the timing and size of analog earthquakes. *Geophysical Research Letters*, 46(3), 1303–1311. https://doi.org/10.1029/2018gl081251

Hulbert, C., Rouet-Leduc, B., Johnson, P. A., Ren, C. X., Rivière, J., Bolton, D. C., & Marone, C. (2019). Similarity of fast and slow earthquakes illuminated by machine learning. *Nature Geoscience*, 12(1), 69–74. https://doi.org/10.1038/s41561-018-0272-8

Jasperson, H., Bolton, D. C., Johnson, P., Guyer, R., Marone, C., & de Hoop, M. V. (2021). Attention network forecasts time-to-failure in laboratory shear experiments. *Journal of Geophysical Research: Solid Earth*, 126(11), e2021JB022195. https://doi.org/10.1029/2021jb022195

Johnson, P. A., Rouet-Leduc, B., Pyrak-Nolte, L. J., Beroza, G. C., Marone, C. J., Hulbert, C., et al. (2021). Laboratory earthquake forecasting: A machine learning competition. *Proceedings of the National Academy of Sciences*, 118(5), e2011362118. https://doi.org/10.1073/pnas. 2011362118

Kachanov, M., & Sevostianov, I. (2012). Rice's internal variables formalism and its implications for the elastic and conductive properties of cracked materials, and for the attempts to relate strength to stiffness. *Journal of Applied Mechanics*, 79(3), 031002. https://doi.org/10.1115/1. 4005957

Kandula, N., Cordonnier, B., Boller, E., Weiss, J., Dysthe, D. K., & Renard, F. (2019). Dynamics of microscale precursors during brittle compressive failure in Carrara marble. *Journal of Geophysical Research: Solid Earth*, 124(6), 6121–6139. https://doi.org/10.1029/ 2019jb017381

Karimpouli, S., Kwiatek, G., Martínez-Garzón, P., Dresen, G., & Bohnhoff, M. (2024). Unsupervised clustering of catalog-driven features for characterizing temporal evolution of labquake stress. *Geophysical Journal International*, 237(2), 755–771. https://doi.org/10.1029/ jb082i030p04893

Kranz, R., & Scholz, C. (1977). Critical dilatant volume of rocks at the onset of tertiary creep. *Journal of Geophysical Research*, 82(30), 4893– 4898. https://doi.org/10.1029/jb082i030p04893

Laurenti, L., Tinti, E., Galasso, F., Franco, L., & Marone, C. (2022). Deep learning for laboratory earthquake prediction and autoregressive forecasting of fault zone stress. *Earth and Planetary Science Letters*, 598, 117825. https://doi.org/10.1016/j.epsl.2022.117825

Lockner, D., Byerlee, J., Kuksenko, V., Ponomarev, A., & Sidorin, A. (1991). Quasi-static fault growth and shear fracture energy in granite. *Nature*, 350(6313), 39–42. https://doi.org/10.1038/350039a0

Mattéo, L., Manighetti, I., Tarabalka, Y., Gaucel, J.-M., van den Ende, M., Mercier, A., et al. (2021). Automatic fault mapping in remote optical images and topographic data with deep learning. *Journal of Geophysical Research: Solid Earth*, 126(4), e2020JB021269. https://doi.org/10. 1029/2020jb021269

McBeck, J., Cordonnier, B., Cooke, M., Fattaruso, L., & Renard, F. (2023). Deformation evolves from shear to extensile in rocks due to energy optimization. *Communications Earth & Environment*, 4(1), 352. https://doi.org/10.1038/s43247-023-01023-w

McBeck, J. A., Aiken, J. M., Mathiesen, J., Ben-Zion, Y., & Renard, F. (2020). Deformation precursors to catastrophic failure in rocks. *Geophysical Research Letters*, 47(24), e2020GL090255. https://doi.org/10.1029/2020gl090255

McBeck, J. A., Cordonnier, B., & Renard, F. (2021). The influence of spatial resolution and noise on fracture network properties calculated from X-ray microtomography data. *International Journal of Rock Mechanics and Mining Sciences*, 147, 104922. https://doi.org/10.1016/j.ijrmms. 2021.104922

McCormick, C. A., & Rutter, E. H. (2022). An experimental study of the transition from tensile failure to shear failure in Carrara marble and solnhofen limestone: Does "hybrid failure" exist? *Tectonophysics*, 844, 229623. https://doi.org/10.1016/j.tecto.2022.229623

Paterson, M. S., & Wong, T.-F. (2005). *Experimental rock deformation: The brittle field* (Vol. 348). Springer.

Peng, S., & Johnson, A. (1972). Crack growth and faulting in cylindrical specimens of Chelmsford granite. In *International Journal of rock mechanics and mining sciences & geomechanics abstracts* (Vol. 9, pp. 37–86). https://doi.org/10.1016/0148-9062(72)90050-2

Qi, C. R., Su, H., Nieaaner, M., Dai, A., Yan, M., & Guibas, L. J. (2016). Volumetric and multi-view CNNs for object classification on 3D data. In *Proceedings of the IEEE conference on computer vision and pattern recognition* (pp. 5648–5656).

Reches, Z., & Lockner, D. A. (1994). Nucleation and growth of faults in brittle rocks. *Journal of Geophysical Research*, 99(B9), 18159–18173. https://doi.org/10.1029/94jb00115

Renard, F. (2017). Critical evolution of damage towards system size failure in a crystalline rock [Dataset]. *Nird Research Data Archive*. https:// doi.org/10.11582/2017.00025

Renard, F. (2018). Volumetric and shear processes in crystalline rock during the approach to faulting [Dataset]. *Nird Research Data Archive*. https://doi.org/10.11582/2018.00023

Renard, F. (2019). Dynamics of microscale precursors establish brittle compressive failure in Carrara marble as a critical phenomenon [Dataset]. *AGU Publications*. https://doi.org/10.11582/2019.00001
Renard, F. (2021). X-ray tomography data of westerley granite [Dataset]. Sigma2. https://doi.org/10.11582/2021.00002

Renard, F., Cordonnier, B., Dysthe, D. K., Boller, E., Tafforeau, P., & Rack, A. (2016). A deformation rig for synchrotron microtomography studies of geomaterials under conditions down to 10 km depth in the earth. *Journal of Synchrotron Radiation*, 23(4), 1030–1034. https://doi.org/10.1107/s1600577516008730

Renard, F., McBeck, J., Kandula, N., Cordonnier, B., Meakin, P., & Ben-Zion, Y. (2019). Volumetric and shear processes in crystalline rock approaching faulting. *Proceedings of the National Academy of Sciences*, 116(33), 16234–16239. https://doi.org/10.1073/pnas.1902994116

Renard, F., Weiss, J., Mathiesen, J., Ben-Zion, Y., Kandula, N., & Cordonnier, B. (2018). Critical evolution of damage toward system-size failure in crystalline rock. *Journal of Geophysical Research: Solid Earth*, 123(2), 1969–1986. https://doi.org/10.1002/2017jb014964

Roeloffs, E. A. (1988). Hydrologic precursors to earthquakes: A review. *Pure and Applied Geophysics*, 126(2–4), 177–209. https://doi.org/10.1007/bf00878996

Ross, Z. E., Ben-Zion, Y., & Zaliapin, I. (2022). Geometrical properties of seismicity in California. *Geophysical Journal International*, 231(1), 493–504. https://doi.org/10.1093/gji/gga189

Rouet-Leduc, B., Hulbert, C., Lubbers, N., Barros, K., Humphreys, C. J., & Johnson, P. A. (2017). Machine learning predicts laboratory earthquakes. *Geophysical Research Letters*, 44(18), 9276–9282. https://doi.org/10.1002/2017gl074677

Rudnicki, J. W., & Rice, J. (1975). Conditions for the localization of deformation in pressure-sensitive dilatant materials. *Journal of the Mechanics and Physics of Solids*, 23(6), 371–394. https://doi.org/10.1016/0022-5096(75)90001-0

Scholz, C. (1968). Experimental study of the fracturing process in brittle rock. *Journal of Geophysical Research*, 73(4), 1447–1454. https://doi.org/10.1029/jb073i004p01447

Scholz, C. H. (2019). *The mechanics of earthquakes and faulting*. Cambridge University Press.

Scholz, C. H., Sykes, L., & Aggarwal, Y. (1973). Earthquake prediction: A physical basis. *Science*, 181, 4102–4810. https://doi.org/10.1126/science.181.4102.803

Shreedharan, S., Bolton, D. C., Rivière, J., & Marone, C. (2021). Machine learning predicts the timing and shear stress evolution of lab earthquakes using active seismic monitoring of fault zone processes. *Journal of Geophysical Research: Solid Earth*, 126(7), e2020JB021588. https://doi.org/10.1029/2020jb021588

Su, H., Maji, S., Kalogerakis, E., & Learned-Miller, E. (2015). Multi-view convolutional neural networks for 3D shape recognition. In *Proceedings of the IEEE international conference on computer vision* (pp. 945–953).

Wang, K., Johnson, C. W., Bennett, K. C., & Johnson, P. A. (2022). Predicting future laboratory fault friction through deep learning transformer models. *Geophysical Research Letters*, 49(19), e2022GL098233. https://doi.org/10.1029/2022gl098233

Zaliapin, I., & Ben-Zion, Y. (2022). Perspectives on clustering and declustering of earthquakes. *Seismological Research Letters*, 93(1), 386–401. https://doi.org/10.1785/0220210127