Machine Learning Predicts the Timing and Shear Stress Evolution of Lab Earthquakes Using Active Seismic Monitoring of Fault Zone Processes.

Machine learning (ML) techniques have become increasingly important in seismology and earthquake science. Lab-based studies have used acoustic emission data to predict time-to-failure and stress state, and in a few cases, the same approach has been used for field data. However, the underlying physical mechanisms that allow lab earthquake prediction and seismic forecasting remain poorly resolved. Here, we address this knowledge gap by coupling active-source seismic data, which probe asperity-scale processes, with ML methods. We show that elastic waves passing through the lab fault zone contain information that can predict the full spectrum of labquakes from slow slip instabilities to highly aperiodic events. The ML methods utilize systematic changes in P-wave amplitude and velocity to accurately predict the timing and shear stress during labquakes. The ML predictions improve in accuracy closer to fault failure, demonstrating that the predictive power of the ultrasonic signals improves as the fault approaches failure. Our results demonstrate that the relationship between the ultrasonic parameters and fault slip rate, and in turn, the systematically evolving real area of contact and asperity stiffness allow the gradient boosting algorithm to "learn" about the state of the fault and its proximity to failure. Broadly, our results demonstrate the utility of physics-informed ML in forecasting the imminence of fault slip at the laboratory scale, which may have important implications for earthquake mechanics in nature.

Introduction

Machine learning (ML) methods have followed improvements in geophysical techniques, instrumentation, and data availability over the past decade to rapidly emerge as indispensable toolkits for the geophysical community (Bergen et al., 2019; Kong et al., 2019). For instance, significant effort has been devoted to using ML to improve event detection, arrival‐time picking, phase association, and earthquake location (McBrearty et al., 2019; Mousavi et al., 2020; Perol et al., 2018; Ross et al., 2018; Trugman & Ross, 2019; Wu et al., 2018; Yoon et al., 2015; Zhu & Beroza, 2019). Importantly, this has seen a revival in earthquake forecasting research, particularly focused on applying ML techniques to lab data on fault friction (Rouet‐Leduc et al., 2017, 2018) and rock damage (McBeck et al., 2020) to infer failure modes and predictability (Corbi et al., 2019). Notably, recent studies have successfully demonstrated that both cataloged (Lubbers et al., 2018) and continuous (Hulbert et al., 2019; Rouet‐Leduc et al., 2017) acoustic emission (AE) data can be used to infer fault friction and predict the timing, shear stress, and in some cases the size of labquakes. Moreover, ML has been used to demonstrate that slow and fast earthquakes share similar physics (Hulbert et al., 2019). With varying degrees of success, these techniques can predict field observations of volcanic eruption (Ren et al., 2020) and subduction zone fault slip (Corbi et al., 2020; Hulbert et al., 2020; Rouet‐Leduc et al., 2019).

Broadly, ML techniques can be grouped as supervised or unsupervised. The former involves predetermined features that are mapped to labeled data sets in order to construct a regression model that typically involves highly nonlinear functions. Unsupervised learning models are used where such labeled data sets may not be available or the focus is on trying to identify patterns embedded within the data (e.g., Bolton et al., 2019; Olivier et al., 2018; Tan et al., 2006). Especially as applied to laboratory data, supervised ML techniques have relied on the systematic evolution of statistical “features” such as the AE variance or energy (e.g., Rouet‐Leduc et al., 2017, 2018). Indeed, based on the existing works, one could conclude that systematic changes in microseismicity (AE) indicative of fault zone criticality are required to successfully predict labquakes or tectonic fault slip. However, while precursors are routinely documented in laboratory studies of fault failure and brittle fracturing (Brace et al., 1966; Hedayat et al., 2014; Kaproth & Marone 2013; Sammonds et al., 1992; Scholz, 1968; Scuderi et al., 2016; Shreedharan et al., 2020), observations of systematic precursors or foreshocks prior to earthquakes in nature are not routinely documented (Bakun et al., 2005; Main et al., 2012; Niu et al., 2008). Additionally, even when precursory slip is present, it may be masked by changes in elastic properties of the wallrock (Chiarabba et al., 2020; Shreedharan et al., 2021). Moreover, the underlying physics of ML prediction of lab earthquakes is poorly understood. One recent study by Bolton et al. (2020) demonstrated that the precursory increase in AE energy prior to fault failure is likely linked to preseismic fault slip. Thus, the question arises whether detection of preseismic fault slip and lab foreshocks are necessary conditions for ML‐based prediction of lab earthquakes, and ultimately, whether the same is true for field observations.

Here, we use high‐resolution time‐lapse active seismic monitoring to document the evolution of P‐wave amplitude and velocity during the lab seismic cycle of highly variable and aperiodic labquakes. We find that the gradient boosting the ML algorithm (Friedman, 2001) can accurately predict both the timing and shear stress state of labquakes using a subset of the amplitude and velocity features. Interestingly, our results show that the ML predictive power improves as the fault approaches failure. Our previous works and those of others have established the physical links between systematic variations in ultrasonic attributes (P‐wave amplitudes and velocities), fault zone preslip, and wallrock stiffening throughout the laboratory seismic cycle (Shreedharan et al., 2020, 2021). Further, these ultrasonic attributes can be linked to asperity scale mechanical parameters such as the real area of contact and asperity stiffness (Hedayat et al., 2014; Kendall & Tabor, 1971; Kilgore et al., 2017; Pyrak‐Nolte et al., 1990; Shreedharan et al., 2019). This has a remarkable advantage over the use of AE‐based statistical features whose connections to asperity mechanics are, as yet, poorly understood. Thus, by successfully training our ML algorithm on these active‐source ultrasonic attributes, we assign a physical underpinning for the predictive power of the ML approach.

Methods

Friction Experiments

We performed frictional shear experiments in a servo‐controlled biaxial testing apparatus using a double‐direct shear (DDS) configuration (Karner & Marone, 1998). We sheared two frictional interfaces created by mating three blocks of Westerly granite, with a fine layer (<200 μm thickness and ∼0.25 g/layer by mass) of quartz powder between the interfaces to simulate frictional wear material (Figure 1; also see Shreedharan et al., 2020). The granite surfaces were roughened with #60 grit silicon carbide thus producing a root‐mean‐squared roughness of ∼2 μm. The surface roughness and the dusting of quartz powder (median particle size of 10.5 μm) have comparable dimensions; thus frictional processes include direct contact of the wallrock, wear, and internal deformation in the gouge. The fault normal stress was held constant via a fast‐acting servocontroller. The biaxial testing apparatus is fully servocontrolled with independent hydraulic pistons supplying normal and shear loads. The experiments were instrumented with calibrated strain‐gauge load cells to measure normal and shear stresses, and direct current differential transformers (DCDTs) to measure fault normal and shear displacement. An additional DCDT was attached to the center block close to one of the frictional interfaces, to measure fault slip (Figure 1). All load cells used in this study have a resolution of ±5 N and the DCDTs have a displacement resolution of ±0.1 μm. We acquired mechanical data at 10 kHz and averaged in real‐time to 100 or 1,000 Hz prior to recording. A constant shear displacement rate was prescribed for the longer, central block of the DDS configuration. In all experiments reported here, the samples had a constant nominal frictional contact area of 25 cm2 and were subjected to a 10‐MPa normal stress. The prescribed background shear rate was set at 11 μm/s. All experiments were performed at a nominal room temperature range of 22–24°C and the sample humidity was maintained at 100% to ensure reproducibility. We produced a spectrum of slip modes from slow to fast (Inset to Figure 1; Leeman et al., 2016) by varying the stiffness of the shear loading system (Shreedharan et al., 2020).

Figure 1

Data for one complete experiment showing the shear stress evolution as a function of load point displacement. Unload‐reload cycles at ∼2 and 3 mm promote shear localization and steady state shear. Training and test data come from 18–20 mm to 20–21 mm of shear displacement, respectively. Insets show: Schematic of the experiment setup, zoom of slow, and fast labquakes with definitions of time‐since‐failure (T prev) and time‐to‐failure (T next) for a given event, and a sample elastic wave with the wavelet (blue) used to measure P‐wave travel time and amplitude.

Ultrasonic Monitoring

We conducted continuous seismic monitoring via ultrasonic pulses transmitted through the fault interfaces. We used 500 kHz broadband, P‐polarized lead‐zirconate‐titanate (PZT) crystals (Boston Piezo‐Optics Inc. PZT‐5A 0.5” diameter) and transmitted half‐sine pulses every 1 ms continuously throughout frictional shear. The PZT sensors were epoxied in blind holes within steel platens and positioned adjacent to the granite side blocks of the DDS configuration (See inset to Figure 1). Each pulse was sampled at 25 MHz, thus ensuring high temporal resolution in sampling. We use the first P‐wave arrival to calculate travel‐times and velocities, and take the largest peak‐to‐peak amplitude within the first 5 µs for acoustic transmissivity (Inset to Figure 1). Following Nagata et al. (2014), we report transmissivity, |T|, as

In Equation 1, the square root term accounts for the two frictional interfaces that the ultrasonic pulses traverse (Shreedharan et al., 2019). We calculate the velocity as the ratio of the distance traversed by the wavelet to the travel time through the DDS blocks, after accounting for the time spent in steel and at interfaces. We refer readers to Shreedharan et al. (2021) for a detailed description of our velocity and transmissivity measurements.

ML Model – Training and Testing

We analyze the acoustic and mechanical data using ML to predict the temporal evolution of shear stress and the time remaining prior to failure for multiple laboratory seismic cycles (Figure 2). Because we use a supervised learning approach, our first step involves data preparation and selection of features and labeled data sets. Previous studies have demonstrated that the P‐wave amplitudes and velocities evolve systematically throughout the lab seismic cycle. This, in addition to their direct proxy relationship to asperity deformation mechanics, makes them good candidate features for our model. The labeled data, that is, the “unknown” data sets that we want to predict are the fault shear stress, time since the previous failure event and, importantly, the time remaining until failure (see Figure 1 inset). Figures 2a and 2b show the labeled data and Figures 2c and 2d show the features used in this study. Note that while the shear stress is readily determined as part of regular data collection during the experiment, the time‐since‐failure and time‐to‐failure labels are measured from the shear stress data. The time‐since‐failure for an event is computed from the shear stress minimum of the previous event, which has a time‐since‐failure label of 0 s, to the shear stress maximum for the next event. Conversely, the time‐to‐failure for each stick‐slip cycle is counted from the shear stress minimum at the end of the previous cycle, up to the shear stress maximum corresponding to the current cycle, which represents a time to failure of 0 s. We calculate the time‐since‐failure and time‐to‐failure continuously at each point in the seismic cycle for multiple seismic cycles (Figure 2b). During the coseismic portion of the cycle, that is, the stress drop, the time‐since, and time‐to‐failure are set to zero and these data are not used for the ML regression.

Figure 2

Data labels (a, b) and feature engineering (c)–(f) for the ML models. (a) Shear stress over multiple lab seismic cycles showing a complex range of stress drops and recurrence intervals; (b) Time‐since‐ (dashed) and time‐to‐failure (solid) for the events in panel (a); (c) and (d) raw (green) and smoothed (black) P‐wave amplitude (c) and velocity (d); (e) and (f) relative changes in P‐wave amplitude (e) and velocity (f) normalized against the start of each stick‐slip cycle (black) and an offset vector (blue) of the normalized amplitudes.

As part of feature preparation, we smooth the features using a backward‐looking 10‐point moving average to reduce feature‐side noise in the ML regression (see Figure S1 for the ML regression performance without feature denoising). We then normalize the features for each labquake against the minimum value of the feature at the end of the previous cycle (Figures 2e and 2f). This is done to ensure that long‐term trends in the features due to layer‐thinning, shear localization, wear product formation, and smoothening of the granite surfaces do not overprint on their shorter term evolution during the slip instabilities. Further, this ensures that the features in any slip instability are independent and have no “memory” of the fault's previous state in an earlier cycle. Without employing the normalization procedure, we find that the ML models additionally learn the long‐term evolution of the features, which has a significantly larger amplitude. This results in poor fits to the shorter term evolution of the labeled data sets (Figure S2). As a final feature engineering step, we estimate a time‐advanced version of the feature (amplitude or velocity). This is done by shifting the feature vector forward by five time steps (0.05 s). We use the original feature corresponding to time t and a time‐advanced feature corresponding to time t – 0.05 s in our regression. The temporal evolution of the amplitude and velocity features includes an increasing and comparably decreasing gradient. When regressing these features against a monotonically varying label such as shear stress or time to failure, there ceases to exist a one‐to‐one functional mapping, that is, nonunique solutions. This results in poor performance by regression‐based supervised learning techniques, as noted by previous works (Rouet‐Leduc et al., 2017) because, for a given stick‐slip, two data labels can correspond to the same feature value. We conducted a series of benchmarking tests and found that this “offsetting” procedure solves the problem by informing the ML models whether the feature corresponding to a data label comes from the increasing or decreasing gradient space. Additionally, because we use a time‐advanced offset, the new offset feature is backward‐looking in time, thus eliminating any potential data leakage issues. This method is similar to and was derived from the subwindowing procedure utilized by Hulbert et al. (2019) to solve a similar functional mapping problem in their ML regression. We note here that the offsetting is a mathematical transformation rather than a physical one. Hence, the choice of an n‐point offset is arbitrary and this parameter, in general, can and should be optimized to ensure good fits to the labeled data. We selected a 5‐point offset for this study since it provided excellent fits to the test data without overfitting the training data.

Subsequently, we split our labeled data set and features into training and testing sets (Figure 1). While experimental trials included training and test data from displacement ranges across the experiment, we ultimately focus on the displacement range of 18–22 mm in this manuscript (Figure 1). This specific shear displacement range was selected to ensure that comparisons across experiments with different loading stiffnesses are done at similar shear displacements (and strains). Consistent with previous studies, we report on the results from a 70‐30 split in a contiguous fashion, that is, the initial 70% of the data were designated as the training set and the remaining 30% as the test data set (Hulbert et al., 2019; Rouet‐Leduc et al., 2017, 2018). Note that, unlike previous ML studies on experimental faults, which make use of windowing procedures to estimate features (Bolton et al., 2019; Corbi et al., 2019; Hulbert et al., 2019; Rouet‐Leduc et al., 2017), our mechanical and ultrasonic data are synchronized in time. In other words, for each temporally evolving value of the labels (shear stress, time‐since‐failure, and time‐to‐failure), there exists a corresponding data‐point in the feature‐space (amplitudes and velocities).

We use the gradient boosting ML algorithm based on decision trees (Friedman, 2001; Hulbert et al., 2019), to jointly analyze our acoustic and mechanical data sets. Specifically, we utilize an open‐source implementation of this algorithm named XGBoost (Chen & Guestrin, 2016). Hyperparameter tuning is the first step in implementing this model. Hyperparameters define and determine the parameter‐space of functions that can serve as potential models, and the model performance is significantly sensitive to the assigned hyperparameters. The XGBoost implementation has a suite of hyperparameters which must be optimized prior to training. In this study, we determine the optimal hyperparameters by implementing an Efficient Global Optimization (EGO) function to minimize a misfit function (Jones et al., 1998). The optimal hyperparameters are determined via a five‐fold cross‐validation using the K‐Folds cross‐validator (See the supporting information text for additional information and hyperparameter values). Here, a subset of the training data is modeled using a given set of hyperparameters and validated against the remaining training data. Once the hyperparameter “tuning” step is complete, we train the ML model by performing regressions on the feature‐label relationship. This step involves iteratively developing the structure of the gradient boosted decision trees. For a detailed description of this method, we refer the readers to the XGBoost documentation (xgboost.readthedocs.io) and Hulbert et al. (2019). We also calculate the feature importance (see the supporting information text and Figure S3) as an F‐score for the different ML models and find that they are comparable, with the amplitude feature being deemed marginally more important than the velocity. Finally, the optimal model is tested and its performance is compared against the labels, that is, the true experimental values of shear stress and time to failure. We evaluate the model performance in two ways: (a) a qualitative benchmark against a naïve constant‐recurrence interval model for the time‐to‐failure and (b) a quantitative estimation of the model performance using the coefficient of determination (R 2) and the root mean squared error (RMSE) for all labeled data.

Results

Friction Data and Stick‐Slip Periodicity

Our lab results include the full spectrum of frictional slip modes from slow to fast labquakes. In particular, the experimental data are characterized by a range of aperiodic frictional slip instabilities (Figure 3) at the friction stability boundary (e.g., Gu et al., 1984; Leeman et al., 2016). Figure 3a quantifies this frictional chaos by plotting the interevent times as a cross‐plot of time since the last instability T prev and time to the next instability T next (See inset to Figure 1). Data that fall on the 1:1 line represent perfect periodicity, that is, they are time and slip predictable events (e.g., Shimazaki & Nakata, 1980). The 1:2 and 2:1 lines form an envelope representing period doubling (Veedu et al., 2020). Further, data points are colored by the magnitude of stress drop for these events. Broadly, while our data have periodicities between perfectly periodic and doublets, the higher stress drop events are generally represented by more aperiodic behavior while the smaller events cluster tightly around the 1:1 line.

Figure 3

Complexity and periodicity of the labquake recurrence interval. (a) Cross‐plot of the time elapsed since the previous event and time remaining to the next event colored by labquake stress drop. Fiducial lines for 1:2 and 2:1 represent period doubling and 1:1 represents perfectly periodic events. (b) Evolution of peak and minimum shear stress for multiple events in the training and testing data sets, plotted by event number and thus increasing shear displacement. Note that stress drop increases slightly with shear (c) ratio of time since previous and time to next event plotted by event number. The dashed lines at 0.5 and 2 represent the envelopes for period doubling. A ratio of 1 represents perfectly periodic events.

We observe that stick‐slip instability behavior and frictional chaos evolve with shear displacement, represented by increasing event number (Figure 3). Event numbers are calculated from the ML train‐test catalog (Figure 1), that is, between load point displacements of 18 and 21 mm, which contain 220 events. While the peak stress remains relatively constant throughout shear, the stress minimum during the coseismic phase increases with shear (Figure 3b). In other words, the fault starts off in a quasi‐stable condition and becomes increasingly unstable. This is a common observation related to shear driven reduction in the friction critical slip distance (D ) and increase in the friction rate parameter, (b − a), within the rate‐state frictional framework (e.g., Marone, 1998a). The interevent ratio from Figure 3a is plotted in Figure 3c as a function of event number, which is used here as a proxy for shear displacement. For the first ∼50 events, the slip behavior is relatively periodic. As the stress drop increases with shear, the slip behavior becomes increasingly complex which makes this an ideal data set to challenge ML approaches for prediction.

Co‐Evolution of Friction and Elastic Properties

Fault zone elastic properties and frictional strength co‐evolve in a systematic manner during the lab seismic cycle (Figure 4). The fault zone elastic wave amplitude and velocity co‐evolve with shear stress and fault slip rate during both slow and fast slip lab earthquakes (Figure 4a). The corresponding fault displacement, obtained from an onboard displacement sensor (Inset to Figure 1) and slip velocity, measured from the time derivative of this fault slip, also increase systematically as the fault approaches failure (Figure 4b). For our loading rate of 11 μm/s, the fault experiences peak slip rates of ∼300 μm/s, and interseismic locking rates of under 0.1 μm/s for the largest, fastest events. The corresponding transmissivity and P‐wave velocities are shown in Figures 4c and 4d. Note the strong inverse correlation between the fault slip rate and the ultrasonic attributes, which is consistent with observations from previous ultrasonic studies of laboratory stick‐slips (Kilgore et al., 2017; Nagata et al., 2014; Shreedharan et al., 2020, 2021).

Figure 4

Evolution of the mechanical and ultrasonic data over multiple laboratory seismic cycles. Temporal evolution of (a) shear stress (b) slip (black) and slip rate (blue) estimated from the onboard slip sensor. Note that the fault is nearly locked for a big fraction of the seismic cycle (c) transmissivity and (d) P‐wave velocity. Note that both transmitted wave amplitude (c) and velocity show clear precursory changes prior to failure.

Machine Learning Models of Shear Stress, Time‐Since, and Time‐to‐Failure

We document the training, testing and performance metrics of the XGBoost models for shear stress, time‐since‐failure, and time‐to‐failure labels in Figures 5, 6, 7. We train and test two ML models: one for transmitted wave amplitude and one for velocity. We do this to test whether the combination of both physics‐based features (amplitude/velocity) is essential to the performance of our ML models or whether either feature is sufficient for successful model prediction. We do not use a separate validation set; rather we perform a five‐fold cross‐validation on the training data sets.

Figure 5

Shear stress prediction using machine learning (ML). (a)–(c) Prediction using amplitudes as features showing (a) training (b) testing data set, and (c) model performance expressed as a cross‐plot of experimental shear stress versus ML model results with 1:1 line indicating perfect model accuracy. (d)–(f) Prediction using velocities as features showing (d) training (e) testing data set, and (c) model performance. Detailed comparison for a representative seismic cycle showing data (black) and model (blue) trained on amplitude (g) and velocity (h). In both panels, the lower plot (gray) shows the root mean squared error (RMSE) misfit.

Figure 6

Time‐since‐failure prediction. (a)–(c) Prediction using amplitudes as features showing (a) training (b) testing data set, and (c) model performance expressed as a cross‐plot of experimental shear stress versus machine learning (ML) model results with 1:1 line indicating perfect model accuracy. (d)–(f) Prediction using velocities as features showing (d) training (e) testing data set, and (c) model performance. Detailed comparison for a representative seismic cycle showing data (black) and model (blue) trained on amplitude (g) and velocity (h). In both panels, the lower plot (gray) shows the root mean squared error (RMSE) misfit.

Figure 7

Time‐to‐failure prediction. (a)–(c) Prediction using amplitudes as features showing (a) training (b) testing data set, and (c) model performance expressed as a cross‐plot of experimental shear stress versus machine learning (ML) model results with 1:1 line indicating perfect model accuracy. (d)–(f) Prediction using velocities as features showing (d) training (e) testing data set, and (c) model performance Detailed comparison for a representative seismic cycle showing data (black) and model (blue) trained on amplitude (g) and velocity (h). In both panels, the lower plot (gray) shows the root mean squared error (RMSE) misfit.

We report the ML results of shear stress prediction using transmissivity (and its offset) as the sole feature (Figures 5a and 5b) and velocity (and its offset) as the sole feature (Figures 5d and 5e). Here, we quantify the model performance using a standard R 2 metric (Figure 5). Both models (trained on amplitudes and on velocities) have higher R 2 metrics for training compared to the test set (Figure 5), which is normally expected since the models are bound to perform better on data sets they have previously “seen.” Regardless, the test set has a reasonable model performance of R 2 = 0.80 for both amplitude and velocity‐based ML models. Cross‐plots of ML model estimates and experimentally measured values of shear stress show the model performance (Figures 5c and 5f), with the solid (1:1) line representing perfect predictions. We observe that the models tend to deviate from the experimental data early in the seismic cycle, that is, the shear stress minima, and close to the middle of the cycle, approximately where the amplitudes and velocity features switch from an increasing to decreasing gradient. However, the ML model performance improves as time to failure decreases (Figures 5g and 5h). For both amplitude and velocity‐based ML models the RMSE calculated over a 10‐point moving window shows somewhat poor performance early in the lab seismic cycle (i.e., at ∼3,555.4 s) but improved performance later in the cycle, that is, between 3,556 and 3,557 s (Figures 5g and 5h).

Figure 6 shows a snapshot of the results from a training and testing exercise on the time‐since‐failure data label. The ML performance using amplitudes and velocities as features is shown in Figures 6a–6c and  6e–6f, respectively. The ML models perform better during the early parts of the seismic cycle (Figure 6c and 6f). In other words, the time‐since‐failure label is best predicted by the early portions of the amplitude and velocity evolution, immediately following a seismic event. Zooms of representative time‐since‐failure evolution over one seismic cycle, show the corresponding ML model and 10‐point windowed RMSE evolution (Figures 6g and 6h). These plots show the superior fits to the experimental data early in the seismic cycle. Specifically, the RMSE is lowest (or nearly 0) in the initial ∼1 s following a stick‐slip event, as the fault relocks (Figure 4b), and then continually increases until the next stick‐slip event.

ML models focused on the time‐to‐failure label are the most pertinent for earthquake forecasting. Figures 7a–7c show the results from a ML model created using amplitudes as the primary feature, and Figures 7d–7f show the results form a model using velocities. Generally, we document poor performance early in the cycle (see Figures 7c and 7f between 2 and 4 s) when the time‐to‐failure is highest. The model fits improve as the fault approaches failure. Figures 7g and 7h show representative stick‐slip cycles with their associated RMSE for the amplitude‐ and velocity‐trained ML models. Again, we document a reducing RMSE as the fault approaches failure, that is, as the time‐to‐failure approaches zero, particularly in the final ∼0.8 s or final ∼22% of the interseismic period.

Benchmarking Performance Against a Naïve Model

As a final performance evaluation, we benchmark our ML results against a simple model of recurrence times constructed from an averaged recurrence interval for all events in a given experiment (Figure 8). Over the course of the experiment, the (∼200 μm) gouge layer is constantly thinning. Thus, what starts off as a complex interaction between rough granite surfaces and gouge, slowly evolves into an interaction almost exclusively between the granite surfaces as the gouge is slowly sheared out of the system. The average recurrence interval over the entire experiment represents the recurrence of an “average” gouge maturity. In other words, because it assumes a statistical average of all recurrence intervals, this model makes no implicit assumptions about and has no knowledge of the experiment or fault behavior (Figure 8). In our experiments, we see more complex and chaotic periodicity early in the experiments, where there is still significant interaction between the granite and gouge, and a more steady‐state periodicity later as the gouge is removed via shear (Figure 1). With an R 2 = 0.49, the performance of this naïve model is relatively unreliable. When compared with our ML models (Figures 5, 6, 7) which have R 2 metrics of 0.8–0.9 over the entire slip cycle, it is clear that the XGBoost models deliver superior performance.

Figure 8

Naïve model for benchmarking performance. Time to failure for the experimental data (black) and a model (green) based on constant recurrence interval. The low R 2 value shows that the naïve model predicts earthquake failure times poorly.

Discussion

A Physical Basis for the Evolution of Fault Zone Elastic Properties in the Seismic Cycle

The evolution of elastic wave properties around fault zones has been extensively studied in the laboratory (e.g., Passelègue et al., 2018; Paterson & Wong, 2005; Stanchits et al., 2003) and to a smaller degree, on crustal faults (e.g., Brenguier et al., 2008; Chiarabba et al., 2020; Malagnini et al., 2019; Niu et al., 2008). Specifically, systematic variations in the P‐wave velocity field prior to fault failure in the laboratory over multiple slow and fast cycles have been documented (Kaproth & Marone, 2013; Tinti et al., 2016) and attributed to preseismic creep (Scuderi et al., 2016). Similarly, variations in P‐wave amplitudes (or transmissivity or transmission coefficient) have been documented as arising from preseismic creep (Hedayat et al., 2014; Shreedharan et al., 2020). In crustal faults, Malagnini et al. (2019) observed a preseismic attenuation signal in the 20–40 Hz frequency range leading up to the 2004 M6 Parkfield earthquake, which they attributed to fluctuations in the fault zone crack density. However, the ultrasonic amplitudes can be connected to the microscopic asperity stiffness and during stable sliding, the real area of contact (Kendall & Tabor, 1971; Kilgore et al., 2017; Shreedharan et al., 2019). Formalized mathematically, when the ultrasonic pulse wavelength is significantly larger than the fault zone width as is the case here, the specific interface stiffness, k , is related to the transmissivity, |T|, as (Kilgore et al., 2017; Pyrak‐Nolte et al., 1990)

Here, is the density of the surrounding medium (Westerly granite) and is the P‐wave velocity through this medium. More recently, Shreedharan et al. (2021) demonstrated that while the ultrasonic amplitudes track fault creep, seismic velocity contains information about fault creep as well as shear stiffening of the wallrock. In other words, at the onset of microscopic preseismic slip prior to macroscopic failure, the amplitudes start to reduce with an inverse log‐linear relationship with fault (pre‐)slip rate. On the other hand, the ultrasonic velocities increase in the wallrock immediately surrounding the fault due to elastic stiffening, and decrease within the fault due to preslip. This “duality” of information contained in the seismic velocities has been documented in crustal faults as well (Chiarabba et al., 2020).

Based on Equation 2, the strong inverse correlation between wave amplitude and the fault zone slip rate (Figure 4) can be interpreted as the interseismic and coseismic evolution of fault zone asperity stiffness. More specifically, immediately following a stick‐slip event, the fault locks up and heals interseismically (J. H. Dieterich, 1972; Kaproth & Marone, 2014; Marone, 1998b; McLaskey et al., 2012). During this period, the increasing amplitudes and velocities can be interpreted as an increase in asperity stiffness due to reduced slip rate. Similarly, the reduction of precursory amplitude and velocity at and after the onset of preslip can be interpreted as a reduction in asperity stiffness (or, perhaps, asperity destruction) due to welded contact junctions being broken by fault slip. While not a prominent feature of this data set (Figure 4), the amplitudes and velocities are distinctly out of phase due to the additional information pertaining to shear stiffening contained in the velocity (See Figure 2 of Shreedharan et al., 2021). Additionally, the feature importance (F‐score; Supplementary Figure S3) values indicate that the amplitude may be marginally more important than the velocity. This is likely because the ultrasonic amplitudes are stronger markers of fault slip rate (Shreedharan et al., 2020) but the interseismic amplitude and velocity variations are qualitatively similar (Shreedharan et al., 2021). Moreover, this is consistent with the view that preslip is necessary for the production of preseismic AEs (Bolton et al., 2020) which have been used in ML‐based prediction efforts for laboratory instabilities (e.g., Hulbert et al., 2019). In a physical model, this may be translated as an interseismic and preseismic reduction in wallrock crack density due to increased deviatoric stresses experienced by the wallrock during elastic strain energy build‐up. Thus, by training our regression‐based models on transmissivity and velocity, we are, by proxy, training our models on the systematic temporal evolution of the interface stiffness and bulk shear stiffening over multiple seismic cycles.

Connecting ML Model Response to the Physics of Frictional Sliding

The supervised ML approach used in this study “learns” from a regression between the features and labels; thus, no temporal information about a prior state of the fault is explicitly transferred from the experimental data to the ML model. At any point in time, the coefficient of friction, μ, of the fault can be expressed within the rate and state frictional framework as a function of the fault slip velocity, V, and a state evolution term, , as

Here, a and b are the rate‐state constants and the subscript “0” denotes a reference variable. The frictional state term, at steady state is generally thought of as an average asperity lifetime and is described as the ratio of a characteristic slip distance, D , and the fault slip rate, V. The term can only be determined via an iterative inversion procedure by solving Equation 3 simultaneously with a state‐evolution law and an elastic coupling equation. However, Equation 3 can be rearranged (Nagata et al., 2012) to redefine the frictional state as , where

Phase plane cross‐plots between the mechanical attributes of the fault (Figures 9a and 9b), that is, between the friction coefficient, slip rate, and frictional state (estimated using Equation 4 with a = 0.0091), show the inter‐relationships between the rate and state parameters described in Equation 3 (Marone, 1998a). Similarly, cross‐plots of the friction coefficient and the ultrasonic attributes (Figures 9c and 9d) are the relationships studied by the XGBoost models to optimize regressions for shear stress prediction (Figure 5). The phase plane plots in Figure 9 contain hysteretic loops for ∼10 slow and fast seismic cycles. Note the strong inverse correlation between the slip rate (Figure 9a) and the ultrasonic attributes (Figures 9c and 9d). Likewise, the frictional state evolution (Figure 9b) is more directly correlated with these attributes. Qualitatively, our data illustrate the relationships between ultrasonic amplitudes and state suggested by Nagata et al. (2014). This implies a mechanical relationship between the physically determined asperity stiffness and the more empirical frictional state, thought to be an average asperity or contact junction lifetime. However, because our stick‐slip cycles represent the nonsteady state, unstable slip behavior of the fault, no quantitative relationships between these parameters can be derived (Kame et al., 2014; Shreedharan et al., 2019). We also note that, due to the additional gouge in our experimental configuration, some of the shear sliding and instabilities may be accommodated by force chains (Anthony & Marone, 2005; Daniels & Hayman, 2008; Hazeghian & Soroush, 2016; McBeck et al., 2019), particularly earlier in the experiment, that is, at lower shear displacements. Because the granite roughness is comparable to the grain size of the gouge, our shear configuration likely allows granite‐granite shear, granite‐gouge shear, and internal deformation and force chain formation within the gouge. As the experiment proceeds, the gouge layer experiences thinning as the quartz dusting is expelled from the system. Thus, the experiment likely evolves from a force‐chain dominated system early in the experiment to a system that increasingly accommodates real‐area‐of‐contact deformation during the later stages where we perform our ML exercises.

Figure 9

Phase plane cross‐plots of (a) friction and slip velocity expressed as a ratio of the slip velocity to loading rate on a logarithmic scale (b) friction and frictional state expressed on a logarithmic scale (c) friction and P‐wave amplitudes (d) P‐wave velocity. The interseismic locking (green), preseismic creep (blue), and coseismic slip (orange) phases are annotated in each plot.

Here, we consider the potential relationships between the ML model results (Figures 5, 6, 7) and the physical basis for the variations in ultrasonic attributes throughout the seismic cycle (Figure 9e). For a multicontact interface composed of numerous contact junctions, the interseismic healing phase is marked by increasing contact junction size and number of contacts. During this period, the ultrasonic amplitudes and velocities increase (Ryan et al., 2018). Subsequently, as the fault creeps prior to failure, some of these contact junctions are destroyed and shrink in size or cease to exist. This is marked by a reduction in the transmissivity and fault zone velocity (Shreedharan et al., 2021). Finally, during the coseismic slip phase, a number of contact junctions are broken and the fault slips, releasing the stored strain energy in these asperities.

Recall that, for shear stress prediction (Figure 5) and time‐to‐failure prediction (Figure 7), the models perform somewhat poorly during the early portions of the seismic cycle, whereas the model performance is remarkably accurate as the fault approaches failure. In the case of time‐since‐failure prediction (Figure 6), this trend is reversed and the ML model predictions are excellent during the early portions of the seismic cycle, immediately following a stick‐slip event. This is consistent with the observations of Lubbers et al. (2018) who documented that the cataloged AEs early in the seismic cycle were better predictors of time‐since‐failure than time‐to‐failure. Surprisingly, our observations are in contrast with the model of Hulbert et al. (2019). They document an inverse relationship between the duration of the next slip event and the acoustic energy early in the slip cycle of the current event. We report data at the same normal stress, and in turn, the same fractional asperity contact area (Shreedharan et al., 2019) for different stick‐slip magnitudes whereas Hulbert et al. (2019) observe this inverse relationship on data collected at different normal stresses (and thus, different real contact areas). Because stick‐slip magnitudes (Leeman et al., 2016) and AE amplitudes (and energy) scale with normal stress (Bolton et al., 2020; Rivière et al., 2018), the different trends documented here and by Hulbert et al. (2019) could, in part, be explained by normal stress dependence (or a lack thereof) in the experimental designs. Moreover, the relationship between dynamics of AEs and active seismic data are unknown, which makes it harder to compare the results directly. However, we note that Hulbert et al. (2019) document better model fits to their shear stress data closer to failure rather than early in the seismic cycle. Thus, our results are qualitatively consistent in this regard. Based on our observations of the accurate model fits early in the cycle for time‐since‐failure and later in the cycle for time‐to‐failure, we posit that the early, interseismic ultrasonic attributes have significant predictive power and may contain information about the past state of the fault (Interseismic healing in Figure 9e). This is qualitatively similar to the idea that aftershock duration may be linked to mainshock size (J. Dieterich, 1994; Lubbers et al., 2018) although no evidence for this assertion appears to exist in crustal faults (Ziv, 2006). Likewise, the ultrasonic attributes after the onset of preslip likely contain predictive power about the future state of the fault (Preseismic creep phase in Figure 9e), that is, timing and size of the imminent failure (Rouet‐Leduc et al., 2017; Shreedharan et al., 2020; Shokouhi et al., 2021).

For each of the ML feature sets we studied, the model performance suffers during the middle of the seismic cycle, as the fault begins to unlock. We propose that this occurs because the frictional asperities experience a competition between healing, which strengthens contact junctions, and preslip, which rejuvenates contacts. These competing mechanisms translate into nearly zero rate‐of‐change for the features (Figure 4) at the onset of preslip, which in turn, results in poor regression fits because the monotonically varying data labels are associated with nearly constant features. For instance, in the limited number of stick‐slip instabilities plotted in Figures 9c and 9d, notice that the P‐wave amplitudes and velocities are relatively constant at ∼71 a.u. and 5,719 m/s, respectively, for friction coefficients in the range of 0.575–0.59, that is, during the transition from interseismic healing to preseismic creep.

Taken together, our ML model predictions (Figures 5, 6, 7) and the physical mechanisms behind the evolution of the ultrasonic parameters (Figures 4 and 9) indicate that even a relatively simple regression, incorporating no explicit temporal information, can provide insights into the mechanics of fault stability and earthquake nucleation, at least on the laboratory scale. More specifically, because the ultrasonic (amplitude/velocity) evolution is physically related to healing and preslip‐driven interfacial asperity stiffness variations, and shear stiffening of the wallrock, the ML model is able to infer the state of the fault and proximity to failure from these microphysical variations.

These results have important implications for faulting and earthquake forecasting on the crustal scale, particularly near shallow subduction zones (Reasenberg, 1999), and in regions where elastodynamic foreshocks (Dodge et al., 1996; Ellsworth & Bulut, 2018) and slow aseismic creep fronts (Melbourne & Webb, 2002) have been detected prior to a mainshock. Preseismic crustal velocity anomalies have been observed prior to a limited number of earthquakes (e.g., Chiarabba et al., 2020; Niu et al., 2008). Thus, these ML methods could potentially be applied to such regions and, particularly around repeating earthquakes or shallow slow earthquakes, if velocity or attenuation trends are available over multiple seismic cycles.

Conclusions

We study the feasibility of predicting the timing and shear stress of laboratory earthquakes using high‐resolution measurements of transmitted wave amplitude and velocity. Our data provide the first test of using active source seismic data to predict labquakes. We find that supervised ML is capable of predicting the timing and shear stress state of labquakes with reasonable accuracy. Moreover, our results indicate that postseismic increases in the ultrasonic amplitudes and velocities, often associated with fault and frictional healing, may contain memory of the past state of the fault. Importantly, our predictions of fault time‐to‐failure improve in accuracy prior to failure, indicating that fault preslip, which reduces ultrasonic amplitude and velocity, has significant predictive power in the context of imminent failure. Finally, the physical underpinning of the systematic changes in ultrasonic attributes is grounded in the deformation mechanics and the evolution of stiffness of microscopic load‐bearing asperities. Hence, we are able to assign a physical model for the results of our ML models. Overall, our study demonstrates the utility of ML techniques in the study of fault mechanics at the laboratory scale and serves to motivate future pursuits in the quest to improve earthquake forecasting and hazard preparedness in crustal faults.

Supporting Information S1

Click here for additional data file.