Predicting the Timing of Catastrophic Failure During Triaxial Compression: Insights from Discrete Element Method Simulations

Examining the deformation of rocks during triaxial compression may provide insights into the precursory phase that leads to large earthquakes by revealing the components of the deformation field that evolve with predictable, systematic behavior preceding catastrophic failure. Here, we build three-dimensional discrete element method simulations of the triaxial compression of rock cores that include a variety of fault geometries in order to identify the components of the deformation field, including the velocity and strain components, that enable machine learning models to predict the timing of macroscopic failure. The results suggest that the velocity field provides more valuable information about the timing of macroscopic failure than the strain field, and in particular, the velocity component parallel to the maximum compression direction. The models also strongly depend on the second invariant of the strain deviator tensor, J2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${J}_{2}$$\end{document}, indicative of the shear strain. The importance of J2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$${J}_{2}$$\end{document} on the model predictions increases with confining stress, consistent with laboratory observations that show a transition from tensile- to shear-dominated deformation with increasing confining stress. In contrast to expectations from previous machine learning analyses, none of the models strongly depend on components of the strain tensor indicative of dilation, such as the first invariant of the strain tensor. This difference may arise because the simulations host less dilative strain than the experimental rocks.


Introduction
Identifying the factors that control macroscopic failure and how they evolve toward failure are critical goals in geoscience and engineering. Determining which strain components, for example, develop with a systematic evolution towards macroscopic failure may enable analysts to recognize the precursory signals that indicate approaching macroscopic failure in the crust and engineered structures. Precursory activity has been detected preceding some large earthquakes (e.g., Whitcomb et al., 1973). However, for many large earthquakes, including earthquakes that occurred on faults that were heavily instrumented, analyses detected scant evidence of precursors (e.g., Bakun et al., 2005). Variations in seismic wave properties, including the velocity and attenuation, are one of the identified precursory changes in geophysical properties documented in laboratory experiments and some natural earthquakes (e.g., Main et al., 2012). Previous work attributes this variation to the dilation of the crust via fracture propagation and/or opening, preslip on existing fault surfaces, increasing normal stress in the wall rock surrounding the fault, and the movement of fluids in the crust (Main et al., 2012;Scholz et al., 1973). Using the evolution of seismic wave properties, previous machine learning analyses have successfully predicted the timing of slip within precut blocks in the laboratory (Shreedharan et al., 2021a).
The observed reduction in the P-wave velocity suggests that dilation occurs preceding fault slip in the laboratory and the crust, and thus that the dilation of the crust may indicate the proximity of macroscopic failure. However, the evolutions of these seismic properties cannot reveal whether the shear strain accumulates with a recognizable evolution that indicates approaching macroscopic failure. Thus, tracking the seismic properties cannot provide insight into which strain component provides more systematic information, and thus potentially more predictive ability.
To directly compare the predictive ability of the dilation and shear strain, in a previous analysis we used the local incremental strain fields calculated from in situ X-ray tomograms acquired during triaxial compression to predict the timing of macroscopic failure, when the rock experiences a sudden, catastrophic stress drop after achieving its maximum differential stress (McBeck et al., 2020a). This work found that the machine learning models depend primarily on the local dilation of the rock to predict the timing of failure, rather than a quantity that measures the local shear strain components, the curl of the displacement field. This result is consistent with observations from the laboratory and the crust of temporary variations in the P-wave properties prior to fault slip and large earthquakes (e.g., Aggarwal et al., 1973;Chiarabba et al., 2020;Niu et al., 2008;Pio Lucente et al., 2010;Scuderi et al., 2016;Whitcomb et al., 1973). However, we performed this machine learning analysis on experiments with intact rocks that lack a preexisting fault, unlike the laboratory and crustal observations. Because faults can heal and thus strengthen in the interseismic period (e.g., Tenthorey & Cox, 2006), the nucleation process leading to large earthquakes may include fracture propagation and opening within relatively intact rocks, similar to the triaxial compression experiments of intact rocks. However, recent analyses highlight that the variation in seismic properties observed in the crust arises from a combination of processes within the host rock and on the preexisting fault surface that are integrated into the overall signal (Shreedharan et al., 2021b). In laboratory experiments with two sliding blocks, fault preslip tends to reduce the P-wave velocity, while increasing normal stress in the wall rock can increase the P-wave velocity, thereby reducing the strength of any observed reduction of the P-wave velocity (Shreedharan et al., 2021b). Thus, systems that include a preexisting healed fault surface may more closely represent the systems in which earthquakes nucleate in the crust than intact rocks. Moreover, the relative importance of the strain components in predicting the timing of macroscopic failure may differ in systems that include and exclude a preexisting fault. The inclusion of a preexisting healed fault with an orientation that maximizes the Coulomb shear stress will likely promote shear deformation rather than dilation.
To constrain the characteristics of the deformation field that evolve in a systematic evolution toward failure, and thus may be used to recognize the precursory phase of approaching catastrophic failure, we build discrete element method simulations of the triaxial compression of nominally intact rocks, and rock cores that include a preexisting, healed fault. We develop machine learning models to predict the timing of macroscopic failure using characteristics of the deformation field, including the components of the three-dimensional velocity and strain fields. This analysis thus enables comparing the predictive power of each velocity and strain component to each other in systems that more closely resemble natural crustal fault systems than the previous machine learning analysis that used intact rocks (McBeck et al., 2020a). Performing this analysis on simulations, rather than rocks deformed in the laboratory, enables testing a wider range of confining stresses and fault geometries than feasible in the laboratory. We systematically vary the confining stress from uniaxial compression conditions to 50 MPa, vary the orientation of the preexisting fault from 15°to 60°from the maximum compression direction, and vary the fault length from system-spanning to shorter lengths. We test the influence of confining stress on the predictive power of the deformation components because laboratory observations document how confining stress controls the fracture geometry, and proportion of shear deformation relative to dilative or contractive deformation (e.g., Paterson & Wong, 2005). Similarly, laboratory experiments show that the orientation of precut faults can control how the rock fails, and in particular, whether the deformation is concentrated along the preexisting fault surface, or spread more diffusely in the surrounding rock volume (e.g., Giorgetti et al., 2019). We vary the fault length in the simulations because whether the fault spans the rock core, or only crosses a portion of the system, likely controls the precise geometry of the fracture network (e.g., Bobet & Einstein, 1998).
The results suggest that the components of the velocity fields are more valuable for predicting the timing of macroscopic failure than the components of the strain fields, contrary to expectations. In particular, the deformation parallel to the maximum compression direction provides the most useful information to the machine learning models. Some of the machine learning models also depend on the second invariant of the strain tensor, J 2 . Increasing the confining stress increases the dependence of the models on J 2 . None of the models strongly depends on the components of the strain tensor indicative of dilation, the first invariant of the strain tensor, I 1 , even when the models only have access to the normal strain tensor components and I 1 , in contrast to results from previous machine learning analyses. This work thus examines the influence of fault geometry and confining stress on the predictive power of each component of the deformation field, including both the velocity and strain components.

ESyS-Particle Discrete Element Method Simulations
We use the discrete element method tool ESyS-Particle to build the simulations (Abe et al., 2004). ESyS-Particle approximates the behavior of materials with indestructible spherical particles that may be unbonded or connected to other particles with one or several bonds. The bonds have both shear and normal stiffness. During a simulation, a bond may break if the forces acting on it exceed the Coulomb criterion determined from its prescribed microparameters of cohesion and friction. Unbonded particles experience frictional and elastic forces. Thorough calibration of the microparameters of the bonds and particles is required for the simulations to approximate the mechanical response of natural rock. ESyS-Particle requires a damping parameter that decreases the kinetic energy in order to maintain the stability of the simulations. Thus, deformation and fault growth are quasistatic, similar to the laboratory experiments to which we compare the results of the present analysis (McBeck et al., 2020a).
We simulate the uniaxial and triaxial compression of blocks that include one rough fault embedded in a three-dimensional system. The blocks are 40 mm 9 80 mm 9 10 mm along the x, y, and z-axes respectively (Figs. 1, 2). We chose to restrict the dimension of the block along the z-axis in order to improve processing time. Because the preexisting faults extend from the front to the back of the system, along the z-axis, varying the length of this dimension would not significantly change the local stress field that develops within the system, and thus the resulting fracture network.
We apply a constant normal stress to the four sides of the simulations on the front, back, left and right sides, and then apply constant axial velocity loading. To ensure stability, we apply a linear ramp from zero to the desired normal stress and axial velocity over 40,000 timesteps. We identify the timestep size using the maximum bond stiffness and minimum particle radius, following the Courant condition (e.g., Karimabadi et al., 2005) (e.g., Table S1). The normal stress increases in the first 40,000 timesteps, and then the axial velocity loading increases from 50,000 to 90,000 timesteps. Thus, the simulation achieves the desired confining stress before the onset of the axial velocity loading. To apply the axial velocity loading, we move the top and bottom boundaries of the simulations toward each other at 0.125 mm/s. We apply 0, 10, and 50 MPa of confining stress.
To build the preexisting fault, the simulations have three groups of bonds with different microparameters (cohesion and stiffness) that represent the principal slip zone, the damaged rock surrounding this zone, and the more intact host rock (Fig. 1). We refer to these bonds, and the types of fractures that develop within these materials via bond breakage, as on-fault, damage zone, and off-fault, respectively. The simulations include one fault that may be three different lengths, L: 15 mm long, 20 mm long, or system-spanning, in which the fault extends across the x-axis of the simulation (Fig. 2). These faults have roughness parameters of natural faults, with Hurst exponents parallel and perpendicular to the down-dip direction of 0.6, and an amplitude of one millimeter (Candela et al., 2012). Following the approach of McBeck et al. (2021b), we build the faults by labelling particles with surfaces that lie within one millimeter of a mesh surface that represents the fault as on-fault, and then labelling all the other particles as off-fault. Because the maximum Vol. 179, (2022) Predicting the timing of catastrophic failure during triaxial compression 3627 particle radius is one millimeter, using this distances ensures that the resulting fault surface does not contain gaps. The particle labels then determine the bond labels, which indicate whether the bond represents on-fault, damage zone, or off-fault material. Fault zones consist of a principal slip zone and surrounding damage zone embedded within a more intact host rock (e.g., Faulkner et al., 2010). Based on the relationship between fracture density and mechanical properties (e.g., Paterson & Wong, 2005), these materials are expected to increase in strength and stiffness with increasing distance from the fault core. This range in material parameters is also consistent with seismic imaging that reveals low velocity zones around faults, and increasing velocities in the surrounding host rocks (e.g., Allam et al., 2014). We therefore prescribe three different sets of cohesion, c, and normal stiffness, k n , to the on-fault, damage zone, and off-fault bonds. We vary the values of these microparameters to achieve the desired mechanical properties, including the uniaxial compressive strength, UCS, and Young's modulus, E.
We first find the set of microparameters that produces the UCS and E of intact granite, and then prescribe these microparameters to the off-fault bonds. To find these parameters, we build Geometries with 15 mm long faults and the four tested fault orientations: a particles, b bonds. The simulations are 40 mm 9 80 mm 9 10 mm along the x, y, and z-axes, respectively. a The simulations consist of two types of particles that represent the relatively intact host rock (gray), and the weaker preexisting fault material (black). b Three types of bonds connect the particles. The bonds represent the on-fault (yellow), damage zone (blue-green), and off-fault material (dark blue). The on-fault bonds connect the weak on-fault particles. The off-fault bonds connect the strong off-fault particles. The damage zone bonds connect pairs of on-fault and off-fault particles. The faults extend from the front to the back of the system along the z-axis. The faults have the roughness parameters (Hurst exponents) of natural faults. To build the faults, we label particles within 1 mm of a surface with the prescribed roughness parameters (Hurst exponents and amplitude) as the weak onfault particles. The remaining particles in the system are the off-fault particles. We vary the orientation of the fault, h, between 15 to 60°from the maximum compression direction, r 1 homogeneous simulations in which all the bonds have the same set of microparameters, and then load them in uniaxial compression ( Figure S1). Following this step of the calibration, we identify the set of microparameters that produces UCS and E of 270 MPa and 40 GPa, respectively, within the ranges Fault geometries with the four tested orientations and a 15 mm long faults, b 20 mm long faults, and c system-spanning faults Vol. 179, (2022) Predicting the timing of catastrophic failure during triaxial compression 3629 of intact granite (e.g. Griffiths et al., 2017). Table S1 lists the full suite of microparameters found through calibration.
The strength and stiffness of the principal slip zone of a fault depend on the rate of healing, rate of crustal loading and damage accumulation, and time elapsed since the previous earthquake (e.g., Lyakhovsky et al., 2001). We model systems with a relatively weak fault, so we find the set of microparameters for the on-fault bonds that produce UCS and E of 10 MPa and 4.5 GPa, respectively. Then, we prescribe microparameters to the damage zone bonds between the endmembers of the on-fault and off-fault zones. This process produces increasing UCS and E from the on-fault to the damage zone, and then to the off-fault volume. These simulations thus assume that the crustal loading that reactivates the preexisting healed fault occurs quickly compared to the healing of the surrounding host rock.

Machine Learning Models
Throughout this analysis, we refer to the machine learning models as ''models'', and the discrete element method simulations as ''simulations'' in order to differentiate between these tools. We build the models to predict the timing of macroscopic failure, and thus the peak axial stress, using the particle velocities and corresponding strain fields. We do not include information about the local stress fields because the stress tensor is difficult to measure in the field with high spatial resolution. We provide the three-dimensional velocity and strain fields from timesteps of the simulations to the models, and build the models to predict the normalized distance to macroscopic failure. We define the normalized time to failure as ðt f À tÞ=t f , where t f is the timestep when macroscopic failure occurs and t is the present timestep of the simulation. Because the applied axial velocity loading of the simulations is constant, the timestep is proportional to the cumulative axial strain. Thus, the time to failure calculated from the timesteps is equal to the time to failure calculated from the cumulative axial strain. Following previous laboratory work ( McBeck et al., 2020a), we consider the time when the system achieves the maximum axial stress as the time of macroscopic failure.
We use XGBoost regression models to perform these predictions (Chen & Guestrin, 2016). We chose XGBoost because of its efficiency and accuracy (e.g., Bühlmann & Yu, 2003;Friedman, 2001). We develop different sets of models to predict each target for each numerical simulation. The combination of four fault orientations, three fault lengths, and three levels of applied confining stress, leads to 36 unique simulations, with one simulation for each combination of parameters.
From the velocity and strain fields, we extract a series of features as for the machine learning models. We only use the timesteps of the simulations that occur preceding macroscopic failure. We calculate the mean of the three velocity components of all the particles within cubes with 1 mm long sides at increments of 1000 timesteps throughout the simulation. Using this grid of velocity vectors, we calculate the nine components of the three-dimensional strain tensor. We then calculate quantities that measure the volumetric strain: the divergence or first invariant of the strain tensor, I 1 , and quantities that measure the shear strain: the second invariant of the strain deviator tensor, J 2 , and the curl. Table S2 lists all the features.
This method of feature extraction produces several hundred measurements of each feature at each simulation timestep, and thus at various times preceding macroscopic failure. Each sample of the dataset provided to the machine learning models is unique in its combination of time and position. We split the training and testing dataset so that none of the data in the training dataset occurs in the testing data. In particular, we split the data so that data from the same timestep of the simulation does not occur in both the training and testing datasets, and the testing dataset includes about 20% of the data. Consequently, we split the training and testing data sets so that the models are trained on different portions of time in the simulation than they are tested. To capture the influence of random variations in how the data are split, for each of the 36 unique simulations we build ten different models with different splits of training and testing data.
We estimate the influence of each feature on the predictions using a widely used metric: Shapely Additive Explanations (SHAP) (Lundberg & Lee,3630 J. McBeck Pure Appl. Geophys. 2017). We compare the mean absolute value of the SHAP value (mean |SHAP|) of each feature, and thus focus on the overall influence of that feature on the model prediction. We perform a grid search over the hyperparameter-space to find the best set of hyperparameters for each model (Lundberg & Lee, 2017).

Results
Here, we first describe the fracture networks that develop in all of the simulations, and how the confining stress and fault geometry control the fracture networks that develop. Then, we follow the evolution of the fracture networks throughout the simulations up to macroscopic failure, when the simulation achieves the maximum axial stress. We then compare the ability of the models to predict the timing of failure using the strain and velocity fields of the simulations. We finally identify the characteristics of the strain and velocity fields that control the model predictions, and examine how they evolve toward failure.

Fracture Network Development
The confining stress, fault orientation and fault length play key roles in the final fracture geometry that develops by the end of the simulation, following the peak axial stress and subsequent stress drop (Fig. 3). Systems with longer faults tend to concentrate deformation within the preexisting healed fault (i.e., on-fault), whereas systems with shorter faults develop more fractures outside of the preexisting fault (i.e., off-fault). For example, when the fault orientation, h, is 30°, the volume of off-fault fractures is larger when the fault length, L, is 15 mm, than when L = 20 mm or is system-spanning, at each level of applied confining stress. This relationship also holds for the simulations with h = 45°at each level of confining stress, r 2 . However, when the faults are less well-oriented for shear failure, with h = 60°and h = 15°, this relationship breaks down. Instead, fractures develop throughout the off-fault region with similar densities in simulations with varying fault length at the same r 2 . The similarity of the relationship between fault orientation and the magnitude of off-fault deformation in these simulations and laboratory experiments (Giorgetti et al., 2019) provides evidence that the simulations closely approximate the behavior of natural rocks.
The fault length appears to exert a stronger control on the resulting fracture geometry when the preexisting faults are less well-oriented for shear failure, with h = 60°and h = 15° (Fig. 3). In simulations with h = 30°-45°, the off-fault fractures propagate approximately coplanar to the preexisting faults when L is not system spanning, thereby producing systems similar to those with systemspanning preexisting faults. However, in simulations with h = 60°and h = 15°, similar behavior only occurs when r 2 = 0. Higher confining stress suppresses coplanar propagation from the tips of the preexisting fractures in these simulations.
Confining stress also controls the prevalence of fault-oblique fractures (Fig. 3). In simulations with no confining stress, the off-fault fractures tend to develop at oblique angles to the preexisting fault. In simulations with higher confining stress, the prevalence of these fractures tends to decrease. Simulations with system-spanning faults, for example, host either one or no fault-oblique fractures when r 2 [ 0.
Simulations that do not include a preexisting fault develop diffuse fracture networks prior to failure, similar to some laboratory observations (e.g., Renard et al., 2019).
Next, we examine the evolving fracture networks in four characteristic simulations leading to macroscopic failure. These simulations include fault geometries with the longest and shortest tested faults, and orientations that are less and more well-oriented for shear failure, h = 15°and h = 30°. The simulations include three phases of macroscopic stress accumulation due to the applied axial strain (Figs. 4, 5, Figure S2). First, the axial stress increases quasilinearly with the axial strain, then the rate of axial stress accumulation slows, and finally the axial stress decreases. Macroscopic failure thus occurs when the simulations achieve the maximum axial stress.
Fracture development includes several key phases preceding macroscopic failure in simulations with shorter faults, L = 15 mm, (Fig. 4), that do not occur when the faults are system-spanning (Fig. 5). In simulations with L = 15 mm, fractures first develop Vol. 179, (2022) Predicting the timing of catastrophic failure during triaxial compression 3631 along the weaker preexisting fault, and then they propagate from the tips of the preexisting fault (Fig. 4). Throughout this process, fractures also develop off-fault with a diffuse distribution. The confining stress and orientation of the preexisting fault control how the fractures propagate from the tips of the preexisting fault. When the preexisting fault is well-oriented for shear failure, with h = 30°, increasing confining stress promotes coplanar propagation from the tips of the fault. In particular, when the system experiences no confining stress, the fractures propagate approximately parallel to the r 1 direction (Fig. 4c). When the system experiences higher confining stress, the off-fault fractures propagate coplanar to the preexisting fault (Fig. 4d). The observed influence of confining stress on the fracture geometry is consistent with laboratory observations and numerical simulations that observed that increasing confining stress promoted coplanar propagation from the tips of preexisting faults rather than out-ofplane propagation (Madden et al., 2017). When the preexisting fault is less well-oriented for shear failure, with h = 15° (Fig. 4a, b), confining stress has a similar influence on the resulting fracture geometry as when h = 30°. Increasing the confining stress reduces the density of fault-oblique off-fault fractures. Varying h also influences the spatial localization of the off-fault fractures that develop at the tips of the preexisting faults. These fractures are more localized when h = 30°than when h = 15°.
When the faults are system spanning (Fig. 5), fracture development concentrates primarily within the preexisting fault (i.e., on-fault), as expected. Similar to the systems with shorter faults, fracture development first initiates along the preexisting fault, and later fractures propagate throughout the off-fault region. Increasing the confining stress and rotating the fault away from the optimal orientation for shear failure increases the density of off-fault fractures.
The varying spatial distribution of off-fault fractures and propagation from the tips of the preexisting faults due to confining stress, fault orientation and fault length suggest that the characteristics of the velocity and strain fields that the machine learning models use to identify the timing of macroscopic failure may differ between the simulations. For example, the orientations relative to the maximum compression directions of the off-fault fractures that propagate from the tips of the preexisting faults suggest that they accommodate dilation in simulations with L = 15 mm, h = 30°, and r 2 = 0, and accommodate shear strain in simulations with higher confining stress, r 2 = 10 MPa (Fig. 3).
Because we develop the models to predict the timing of macroscopic failure, we now compare the peak axial stress achieved in each simulation (Fig. 6). Increasing the fault length decreases the peak axial stress for simulations with faults well-oriented for shear failure, with h=30°and 45°, at all the tested confining stresses, as expected. Under no confining stress, increasing fault length also promotes lower peak axial stress when the faults are less well-oriented for shear failure, h=15°and 60°. However, at higher confining stress r 2 = 10 MPa and 50 MPa, the peak axial stress is similar for faults of varying lengths when h = 15°and 60°. Thus, when the fault is less well-oriented for shear failure, increasing the fault length does not lead to a significant reduction in the macroscopic strength.
The systems fail at the lowest axial stress when h = 30°and 45° (Fig. 6), as expected from theoretical expectations (e.g., Jaeger & Cook, 1979) and laboratory experiments (e.g., Donath, 1961;Giorgetti et al., 2019). Because these orientations maximize Coulomb shear stress, the lower axial stresses at these orientations and the corresponding approximately parabolic distribution of the peak axial stresses relative to h provides additional evidence that the simulations closely approximate the behavior of natural rocks. When the faults are system-spanning, the peak axial stress with h = 30°is significantly lower than with h = 45°, producing a more asymmetric distribution of the peak axial stress relative to h than when the faults are shorter. Similar to the varying fracture geometries, these varying failure stresses suggest that the machine learning models may use different characteristics of the strain and velocity fields to predict the timing of failure.

Predicting the Proximity to Macroscopic Failure
To assess the ability of the models to predict the proximity to macroscopic failure, we examine the R 2 score of the models. High positive R 2 scores ([ 0.7) indicate strong correlations between the predictions and observations. Scores of one indicate a perfect correlation between the predictions of the models and the observations. The models perform exceptionally well, with mean R 2 scores between 0.93 and 0.98 across all the simulations (Fig. 7). The performance of the models developed for simulations with the highest tested confining stresses is slightly lower than the performance of the other models.
To identify the features that control the model predictions, we calculate a cumulative representation of the SHAP values of each feature. Because we develop ten models for each simulation that only differ in how the training and testing datasets are split, we calculate a normalized importance s=maxðsÞ from the mean |SHAP| value of each feature, s. We then find the mean of the normalized values weighted by the R 2 score of the model across all the models for a given simulation, s s ¼ P R 2 ðs=maxðsÞÞ=n, where n is the number of models. We weight the normalized importance by the R 2 score of the model so that more successful models (with higher R 2 ) have a greater influence on the results compared to less successful models. The value s s thus reflects the distribution of the |SHAP| values across all the features of all the machine learning models for a particular simulation (Fig. 8a). We also report the sum of s s across all the simulations (Fig. 8b). To further constrain the features that most strongly control the model results, we identify the features with s s greater than 25% of the maximum s s of each simulation (Fig. 8c).
The cumulative distribution of s s across all the simulations indicates that four features most strongly control the model results: the y-position in the simulation, the velocity components parallel, v y , and perpendicular, v x , to the maximum compression direction, and the second invariant of the strain deviator tensor, J 2 (Fig. 8b). The distribution of s s for the individual simulations, however, indicates some variations between the simulations (Fig. 8a). Examining the characteristics that have s s greater than 25% of the maximum s s for each simulation (Fig. 8c) indicates that all of the models depend strongly on v y , and nearly all of them depend on y. The models that depend strongly on J 2 tend to be those developed for simulations with higher confining stress. In particular, all of the models developed for the simulations with no applied confining stress do not depend on J 2 , except when h = 30°and when the fault is systemspanning or L = 20 mm. The models that depend on v x tend to be those developed for simulations with h\ 60°. When h = 60°and when the simulations do not include a preexisting fault, v x does not strongly control the model predictions. The only exception to this trend is when L = 15 mm and r 2 = 50 MPa.
To better understand why the models depend on the characteristics identified from the SHAP values (Fig. 8), we now examine the evolution of three characteristics identified as highly important, J 2 , v y , and v x , and one that is not identified as important, I 1 (Fig. 9). We track these characteristics at two positions throughout three simulations, in the middle of the simulations and along the preexisting fault (x = 20 mm, y = 40 mm), and at a position in the lower portion of the system that does not cross the preexisting fault (x = 20 mm, y = 20 mm). At each of these positions, we report the mean value of the characteristic across all of the z-positions.
The evolution of v y of the simulations that do not include a fault highlights that these simulations do not experience an extended macroscopic yielding phase (Fig. 9a). Instead, the local deformation is relatively homogeneous until about 10% to failure, when the v y at the middle and lower positions begin to change, increasing at the middle position, and decreasing at the lower position.
These evolutions help indicate why the models depend on the various characteristics. Each of the example simulations experience a relatively systematic evolution of v y toward macroscopic failure, while I 1 fluctuates more rapidly, sometimes switching between negative and positive values, and thus contractive and dilative strain (Fig. 9). The more systematic evolution of v y compared to I 1 allows the models to build reliable relationships between this characteristic and the timing of macroscopic failure that produce strong correlations between the model predictions and the observed timing of failure.
The varying systematic evolution of v x and J 2 also indicates why some models depend on these characteristics while others do not. The velocity component v x fluctuates rapidly from positive to negative values, indicative of movement to the right and left in the simulation under no confining stress with no preexisting fault (Fig. 9a), and in the simulation with a system-spanning fault, h = 60°and r 2 =50 MPa. (Fig. 9c). In both of these simulations, the models do not depend strongly on v x . Similarly, the model developed with the simulation under no confining stress with no preexisting fault does not depend strongly on J 2 , whereas the models developed with the other two example simulations depend on this characteristic (Fig. 8c). Although the simulation without a preexisting fault experiences a relatively systematic evolution of J 2 , similar in its overall evolution to the other two simulations, the maximum value that J 2 achieves in this simulation at the selected sampling points is two orders of magnitude less than the maximum values in the other simulations. Thus, both the value of the characteristic and the consistency of its evolution toward failure appear to control whether the models use a particular characteristic to predict the distance to macroscopic failure.

Predicting Failure with Subsets of the Deformation Components
To further identify the characteristics that control the model predictions, we now examine the model performance when the models only use subsets of the data. The distribution of the SHAP values indicates that models that use the shear strain component J 2 , Vol. 179, (2022) Predicting the timing of catastrophic failure during triaxial compression and/or the velocity components v y and v x , may produce the most accurate models, while models that use the volumetric strain components, such as I 1 , may be the most inaccurate (Fig. 8). To test this idea, we develop models with three different subsets of the data: the components that describe the shear strain, the components that describe the volumetric and normal strain, and the velocity components (Table 1). Consistent with expectations, the models developed with the velocity components perform better than the models developed with only the sjhear strain components and the models developed with the normal strain components (Fig. 10). In particular, the mean R 2 scores of the models developed with the velocity components are 0.12-0.25 higher than the mean R 2 scores of the models developed with the volumetric strain components (Fig. 10e). Models developed from simulations with lower confining stress have the largest difference in the mean R 2 scores. This result occurs because the models developed for simulations with lower confining stress perform better than models developed for simulations with higher confining stress when they only have Evolution of three features that strongly control the model predictions relative to the normalized time to peak axial stress, ðt f À tÞ=t f : v x , v y and J 2 , and one that does not, I 1 , in three example simulations: a when the system has no fault and r 2 = 0, b when the fault is systemspanning, h = 30°, and r 2 = 0, and c when the fault is system-spanning, h = 60°, and r 2 = 50 MPa, at two example locations: in the middle of the simulation at x = 20 mm, y = 20 mm (dark blue), and in the lower portion of the simulation at x = 20 mm and y = 20 mm (light blue). The curves show the mean values of each of the four characteristics at the selected locations throughout the full thickness of the model, i.e., along the z-axis. The sketches in red in b and c show the preexisting faults, the boundaries of the simulations, and the sampling locations Table 1 Velocity and strain characteristics used in the three types of models that use subsets of the data

Subset Characteristics
Velocity components x; y; z; v x , v y ,v z Volumetric and normal strain components x; y; z; e xx , e yy , e zz ,I 1 Shear strain components x; y; z; e xy , e yz , e zx , J 2 , curl 3638 J. McBeck Pure Appl. Geophys.
access to the velocity components. When the models only use the volumetric and normal strain components, the models developed for simulations with higher confining stress perform better than models developed for simulations with lower confining stress. However, contrary to expectations, the models developed with only the shear strain components do not perform better than the models developed with only the volumetric and normal strain components. Instead, the mean R 2 scores of the models developed from the volumetric and normal strain components are 0.01-0.11 higher than the scores of the models developed from the shear strain components (Fig. 10d). The difference in the mean R 2 scores is the lowest when r 2 = 50 MPa. Thus, when the confining stress is higher, the difference in the performance of the models is minor, and sometimes the models developed with the shear strain components perform better than the models developed with the volumetric strain components.
To better understand the key characteristics that control the three subsets of models, we now examine the cumulative distribution of s s across all the simulations (Fig. 11). Similar to the models developed with all of the strain and velocity components, the models developed with only the velocity components primarily depend on v x , v y and the y position (Fig. 11c). This result is consistent with the high s s of these characteristics for the models developed with all of the data (Fig. 8). Similarly, the models developed with the shear strain components primarily depend on the curl of the velocity field, and on J 2 (Fig. 11a), consistent with the distribution of the s s values of the models developed with all of the data. In contrast, the models that only use the normal and volumetric strain components primarily depend on characteristics that were not identified as highly important in the models that use all of the data. These models primarily depend on the normal strain components: e xx , e yy , and e zz . The higher cumulative s s for e yy compared to the other normal strain components suggest why these models perform better or comparably to the models developed with the shear strain components. In particular, the models developed with all of the data and the models developed with the Figure 10 Comparing the performance of models developed with subsets of the strain and velocity components: the mean AE one standard deviation of the R 2 scores of a models developed with the shear strain components, b models developed with the volumetric and normal strain components, and c models developed with the velocity components, and the difference between the mean R 2 score of each simulation for d models developed with the shear strain components and volumetric strain components, and e models developed with the velocity components and the volumetric strain components. Negative values in d indicate that the models developed with the volumetric strain components perform better than the models developed with the shear strain components. Positive values in e indicate that the models developed with the velocity components perform better than the models developed with the volumetric components Vol. 179, (2022) Predicting the timing of catastrophic failure during triaxial compression 3639 velocity components primarily depend on v y . Because e yy is a function of v y , it is not surprising that when we deprive the models of v y , they are able to use e yy to perform the predictions with similar success.

The Importance of the Displacement Field
We develop machine learning models to predict the proximity to macroscopic failure using the components of the strain and velocity tensors. The models perform exceptionally well, with mean R 2 scores ranging from 0.93 to 0.98 across the simulations (Fig. 7). When the models have access to all of the velocity and strain components, the models primarily depend on the y-position, J 2 , v y , and v x (Fig. 8). The varying fracture networks produced in the different simulations (e.g.,Figs. 3,4,5) suggest that different strain and velocity characteristics may help the models predict the timing of macroscopic failure. Indeed, models developed from simulations with lower confining stress do not depend strongly on J 2 , but only on the velocity components. Higher confining stress produces larger magnitudes of J 2 (Fig. 9), which allows the machine learning models to successfully predict the proximity of macroscopic failure using this measure of the shear strain field. The observed increasing dependence of J 2 with increasing confining stress is consistent with experimental observations that show a transition from tensile-dominated deformation to shear-dominated deformation with increasing confining stress (e.g., Paterson & Wong, 2005). In laboratory experiments, for example, during uniaxial compression, the brittle failure of rocks occurs through the propagation and opening of fractures aligned approximately parallel to the direction of r 1 (Tapponnier & Brace, 1976). During triaxial compression, as the applied confining stress increases, the overall orientation of the dominant fracture network can trend more obliquely to the direction of r 1 , closer to the orientation expected for shear failure (e.g., Renard et al., 2019). The orientation of individual fractures, however, within such system-spanning fault networks may remain oriented parallel to the direction of r 1 , rather than along the apparent macroscopic orientation (McBeck et al., 2021a). Despite this difference between the true orientation of individual fractures and the overall fracture network, such fracture networks can accommodate shear strain (McBeck et al., , 2020bRenard et al., 2019).
The varying importance of J 2 , and consistency of the importance of v y and v x (Fig. 8) highlight that the components of the velocity field, and thus the displacement field, provide more reliable information about the proximity of macroscopic failure than the strain field. Moreover, when we develop models from subsets of the strain and velocity fields, the most successful models are those developed with the velocity components, rather than the shear strain components or volumetric and normal strain components (Fig. 10). One may expect that the strain field may provide more valuable information about the timing of macroscopic failure than the displacement field. For example, fault slip and opening can produce sites of higher shear and dilative strain localized on the fault, producing systematic relationships between applied loading and cumulative strain. Because the displacement field within the host rock surrounding a preexisting fault can evolve in unexpected ways due to slip and opening on rough faults, one may expect that the relationship between the displacements produced by fault movement and loading, and the timing of macroscopic failure may be more difficult for machine learning models to recognize than the systematic relationships of cumulative strain and macroscopic failure. Contrary to expectations, the present analysis suggests that the displacement field also provides critical information about the proximity of macroscopic failure, as well as the components of the strain field (i.e., Fig. 12).

The Importance of the Deformation Field Parallel to the Loading Direction
The s s distribution of the characteristics of the models developed with all of the strain and velocity components (Fig. 8), and this distribution for the models developed with subsets of the components (Fig. 11) highlight that the components of the deformation field that describe the displacement or strain parallel to the r 1 direction are critical for successfully predicting the timing of macroscopic failure (i.e., Fig. 12). In particular, all of the models developed with all of the characteristics depend on v y , the velocity component parallel to the axial loading, while only some of the models depend on v x (Fig. 8). When the models are developed with only the normal and volumetric strain components, they primarily depend on the strain component parallel to the r 1 direction, e yy (Fig. 11). Thus, despite the variety of fracture networks that develop with varying confining stresses and preexisting fault geometries, the deformation field parallel to the dominant loading direction provides the most systematic information of the tested velocity and strain components (e.g., Fig. 9). This systematic information enables the machine learning models to make accurate predictions about the timing of macroscopic failure.
Previous machine learning analyses that predicted the proximity of macroscopic failure in triaxial compression experiments did not examine the individual nine components of the strain tensor, and instead only used the curl of the displacement field and the first invariant of the strain tensor, I 1 , to predict the timing of failure (McBeck et al., 2020a).

Figure 12
Sketch of the key conclusions. The models are the most successful, with the highest R 2 score, when they use the velocity components to predict the timing of macroscopic failure, and are more successful when they use the volumetric and normal strain components than when they use the shear strain components. Thus, tracking the velocity field may provide useful information about the timing of macroscopic failure, as well as tracking the strain field. The models that rely only on the velocity components primarily depend on v y , and secondarily on v x . Thus, the velocity component parallel to the maximum compression direction may provide the most useful information about impending failure. The models that rely only on the volumetric and normal strain components depend on the components in the following order: (1) e yy , (2) e xx , (3) e zz . Thus, the normal strain component that is dominated by axial-parallel contraction, e yy , may be more useful for predicting the timing of failure than the other normal strain components, which host both dilation and contraction in these systems. The models that rely only on the shear strain components depend primarily on J 2 , and secondary on the curl. Thus, efforts to recognize the onset of the precursory stage leading to large earthquakes may benefit from tracking J 2 rather than the curl of the displacement field Vol. 179, (2022) Predicting the timing of catastrophic failure during triaxial compression 3641 This analysis found that the models depended primarily on the dilation, the positive I 1 , rather than the curl or the negative I 1 to predict the timing of failure. In contrast, in the present simulations, the models primarily depend on the velocity and strain components that characterize the deformation parallel to the r 1 direction, v y and e yy , and J 2 in models under higher confining stress. Due to the triaxial compression loading conditions in these simulations, preceding macroscopic failure e yy is generally negative throughout loading, indicative of contraction. Thus, when the models have access to all of the strain and velocity components, they do not use any of the volumetric or normal strain components indicative of dilation, and instead rely on y, v x , v y , and J 2 . When the models have access to only the volumetric strain components, they primarily depend on the normal strain component that is dominated by contraction, e yy . In contrast, the models developed from the triaxial compression experiments primarily depend on the dilation, rather than the curl. In the simulations that include a preexisting fault well-oriented for shear failure, at 30°-45°, we expect that shear deformation may dominate the system, and thus that the models should depend on components of the shear strain. However, even for simulations that lack a preexisting fault, the models do not depend on measurements of the dilation, but primarily on J 2 and v y . These models may not rely on the volumetric or dilative strain, I 1 , because the components of the velocity field and J 2 evolve with a more systematic evolution toward failure than I 1 (e.g., Fig. 9). In addition, the numerical simulations accumulate less dilative strain than the intact rocks deformed in the experiments. Indeed, comparing the proportion of the volume of the experiments and simulations that host dilation, where I 1 [0, out of the total volume throughout each experiment and simulation indicate that the laboratory rocks host 40-60% dilation (McBeck et al., 2020b), whereas the simulations host about 30% (Fig. 13). Comparing this proportion in the final stages of the experiments and simulations, from 50 to 100% of the failure strain reveals an even larger difference in the volume proportion between the experiments and simulations. Thus, in the triaxial compression experiments, more and more of the rock volume experiences dilation as the rocks approach macroscopic failure. In contrast, as the simulations approach macroscopic failure, the proportion of the simulation that hosts dilation does not significantly increase. Whereas the magnitude of dilation provides systematic information that the machine learning models can use to predict the timing of macroscopic failure in the laboratory experiments, the dilative strain does not follow such systematic relationships in the numerical simulations, causing the machine learning models to use other components of the velocity and strain tensor to successfully predict failure. The simulations may experience less dilation than the experiments due to the initial porosity caused by the packing of the spherical particles. Whereas these experimental rocks generally have \ 10% porosity (McBeck et al., 2020b), the initial porosity of these simulations is about 20%, although we use the maximum feasible range of particle radii to achieve the lowest initial porosity under the computational constraints. The higher initial porosity of the simulations promotes more contraction, and enables shear deformation to occur without significant dilation. The lower porosity of the laboratory rocks requires that the dilation of rough fracture surfaces occurs before these surfaces can slide. Indeed, dilative strain is often paired with shear strain in these experiments (McBeck et al., 2020b), highlighting the necessity of dilating before or concurrent with shearing.

Conclusions
We develop machine learning models to predict the timing of macroscopic failure in three-dimensional discrete element method simulations of triaxial compression. Comparing the predictive power of components of the velocity and strain fields in these models identify four characteristics that primarily control the model predictions: the y-position in the simulation, J 2 , v y , and v x . With increasing confining stress, the predictive power of J 2 increases, consistent with laboratory observations that show a transition from tensile-dominated to shear-dominated deformation with increasing confining stress (e.g., Paterson & Wong, 2005). In contrast to expectations from previous machine learning analyses (McBeck et al., 2020a;Shreedharan et al., 2021a), observations of varying seismic wave properties prior to some large earthquakes (e.g., Whitcomb et al., 1973), and laboratory observations of decreasing P-wave velocity prior to stick slip events (e.g., Shreedharan et al., 2021b), the volumetric strain component indicative of dilation, I 1 , does not strongly control the model predictions. The models only depend on the normal strain components that quantify both dilation and contraction in these simulations, e xx and e zz , when we force the models to only use the normal and volumetric strain components. However, even these models primarily depend on the normal strain component that is dominated by contraction, e yy . The difference between the present analysis and previous machine learning analysis on the triaxial compression of rocks (McBeck et al., 2020a) may occur because the rocks deformed in the experiments are dominated by dilation, with increasing volume proportions toward failure (McBeck et al., 2020a(McBeck et al., , 2020b, whereas the numerical simulations do not experience similar strain accumulation (Fig. 13). The models developed with only the velocity components perform better than the models developed with only the volumetric and normal strain components, and better than the models developed with only the shear strain components. This result highlights that the velocity field, as well as the strain field, may provide valuable information about the onset of the precursory phase leading to catastrophic failure, such as earthquakes. In particular, the velocity field parallel to the maximum compression direction may provide the most valuable information about the timing of failure, regardless of the orientation of the preexisting fault.