2. Experimental methods and materials

Experiments were conducted using an ambient pressure, cryogenic temperature, servo-hydraulic biaxial friction apparatus (McCarthy and others, Reference McCarthy, Savage, Koczynski and Nielson2016), with modifications to the insulating cryostat and loading procedure to allow measurement of till (Saltiel and others, Reference Saltiel, McCarthy, Creyts and Savage2021). In this double-direct-shear configuration, a central ice block slides against two stationary side blocks, with 10 mm thick layers of pre-compacted and frozen Matanuska till or Barre granite rock on opposite sides of the ice, such that applied horizontal load is resolved as normal stress (horizontal force over the surface area of each side block, 50 mm × 50 mm) and vertical load as shear stress on the two sliding interfaces (vertical force over the surface area of both side blocks, 2 × 50 mm × 50 mm).

As in Saltiel and others (Reference Saltiel, McCarthy, Creyts and Savage2021), we control temperature with Peltier thermoelectric coolers in front and behind the ice block, as well as circulation of chiller fluid through the side blocks where both temperature and flow rate of chiller fluid were actively controlled to reach and sustain the desired temperature. Resistance Temperature Detectors (RTDs) ported directly behind the till or rock monitor the temperature as close to the sliding interfaces as possible.

In this study, all experiments used ~50 kPa of normal stress with displacement control using a load point velocity of 100 μm s^−m (the slope of the black line on the bottom of Fig. 1b), just over 3 km yr^−r, around the surface velocity of the fastest ice streams (Zoet and others, Reference Zoet2020). Since the load point Linear Variable Differential Transformer (LVDT) only has 20 mm of stroke, the load point was stopped halfway through each experiment and then LVDT was reset to complete the rest of the experimental displacement. In this way, every experiment included a hold of about 60 s during which the shear stress relaxed and then reloaded, usually resulting in the largest stress-drop and AE of each experiment (see Fig. 2). This also provides a measurement of the interface's degree of healing. The cryostat and sample assembly are diagramed in Figure 1a, while additional experimental details are described in supplementary text S1.

Figure 1. (a) Schematic of biaxial cryostat with additions of rubber spring, to decrease loading stiffness, AE sensor frozen into central ice bock (pictured within ice in inset on left) and sample displacement measurement. For more details about apparatus see Saltiel and others (Reference Saltiel, McCarthy, Creyts and Savage2021) and supplementary text S1. (b) An example experiment of measured friction drops (in black on top) and stick–slip sample displacement (in red on the bottom) with the steady load point displacement (in black) for reference. Orange arrow shows an example of how friction drop is measured, stress drop is calculated by multiplying friction drop by the constant normal stress of 50 kPa. Instability was induced by apparatus reaching subcritical stiffness (k_cr). (c) An example AE waveform before processing, from a single stress-drop.

Figure 2. Example experiments of the temperature effect on slip stability for (a) rock and (b) till beds. Each experiment begins with stress-drops but, after a hold (described in section 2 above), with increasing temperature the ice starts to slide stably without sudden friction drops or audible stick-slips. The transition to stable sliding occurs ~ 0°C for the rock experiment. In the till experiment, the stability temperature is reached during the hold, but as it is re-cooled stress-drops do not resume until the temperature is below about −2.5°C. Each estimated transition temperature is highlighted with a thin black horizontal line, but temperature is not measured directly at the ice-bed interface, so the temperature at the interface lags that recorded. The lag time (estimated to be ~100 s) is represented by the yellow region right of the measured temperature. It is also apparent that the till experiment has a higher friction and healing rate (as the friction rose more after hold times of similar duration).

We made three additional modifications to the apparatus from Saltiel and others (Reference Saltiel, McCarthy, Creyts and Savage2021). A Linear Variable Inductance Transducer (LVIT) position sensor measures the sample displacement separate from the loading point's preset displacement. This allowed measurement of displacement in each stress-drop ‘slip’ event as well as slip occurring during ‘stuck’ periods and the timing coincident with stress-drops (Fig. 1b). Here we refer only to mechanical or bulk stress-drops, the stress change over a slip event as measured by our vertical load cell, not to be confused with seismologically derived stress-drops. A piece of rubber was inserted into the loading geometry that effectively reduced the stiffness of the apparatus, reaching critical stiffness and allowing stick–slip instability (Zoet and others, Reference Zoet2020). We estimate the effective apparatus stiffness using the mechanical data's reloading slope between stress-drops, relative to the compression of the loading train including rubber, the load point displacement minus sample displacement (Fig. 1b), giving the apparatus stiffness with the rubber to be ~0.1 kPa μm⁻¹, significantly less stiff than was estimated without the rubber ~1 kPa μm⁻¹ (Saltiel and others, Reference Saltiel, McCarthy, Creyts and Savage2021). Additionally, commercial piezoelectric transducers were frozen into the central ice block, facing the right ice-bed interface, to measure AEs. The sensor has a directionality, measuring waves incident on the transducer face, but waves could also arrive after bouncing around the sample. Recorded later arrivals (Video 1 in supplement) likely represent reflections from the other interface. We tried four different transducer types of varying size and frequency sensitivity, settling on Physical Acoustic's Nano 30^TM miniature AE sensor due to its small size and 125–750 kHz response, covering most of our measured frequency content. All AEs analyzed here were recorded with a single Nano-30.

AEs were recorded using a preamplifier and TiePie^TM HS6 differential digital oscilloscope. To ensure we recorded all relevant spectral content in the waveforms, they were recorded at a sample rate of 100 MHz for 1 ms time windows around each event (Fig. 1c). These oscilloscope settings provided the optimal real-time viewing of waveforms as they were being recorded (Video 1 in supplement), but subsequent analysis showed most of the energy was under 1 MHz, and waveforms were down sampled to 10 MHz. Recordings of continuous acoustic signal without applied shear found electrical noise above 3 MHz, so filtering also helped remove persistent noise sources. The oscilloscope was set in rising-limb trigger mode with trigger amplitude set just above the noise level before slip initiates, such that it did not trigger without an audible stress-drop. Since electrical and other sources of noise can vary, this trigger level was adjusted throughout the experiment to maximize the number of captured events and minimize waveforms of purely noise, but some events were missed, and some events triggered by noise or other AE sources were saved.

3. Stick-slip instability at frozen conditions

Both rock and till experiments were undertaken over a range of temperatures, showing the temperature dependence of instability. Analyzing the temperature dependence of AE characteristics is outside of the scope of this article and will be explored in future work. We found stress-drops only at frozen temperatures, < ~0°C for rock and <~−2.5°C for till beds (Fig. 2). It must be noted that temperatures are approximate since they are measured behind the till/rock. This is accurate for stable temperatures, but there is an unknown time lag before the temperature on the sliding interface reaches the recorded temperature. Additionally, when the temperature probe goes above the pressure melting point, ~0°C, the ice will remain at its pressure melting point. The rock or till layer is 10 mm thick and given thermal diffusivities ~1 mm² s⁻¹, we estimated the lag time to be ~100 s, given by the yellow region after the recorded temperature in Figure 2. The observed temperature dependent stability is consistent with rate-weakening friction shown for till beds at ~ −3°C, but not at the pressure melting point (Saltiel and others, Reference Saltiel, McCarthy, Creyts and Savage2021). We estimate (Eqn (1)) the apparatus stiffness with rubber to be the same order of magnitude as the critical stiffness estimated from separate velocity-step experiments under similar conditions ~0.02 kPa μm⁻¹ (Saltiel and others, Reference Saltiel, McCarthy, Creyts and Savage2021). The factor of five difference is consistent with the error inherit in applying estimations of rate-state friction parameters (b – a, D_c) from a single experiment, as well as in our rough estimation of apparatus stiffness. Past studies of ice on rock friction did not find rate-weakening until lower temperatures, <~−18°C for McCarthy and others (Reference McCarthy, Savage and Nettles2017). In that study, experiments above −18°C which exhibited slight rate-strengthening were undertaken at less than half the slip rate, which could affect the rate-dependence (Zoet and others, Reference Zoet2013; Saltiel and others, Reference Saltiel, McCarthy, Creyts and Savage2021) as well as stability more broadly (Schulson and Fortt, Reference Schulson and Fortt2012). It is also possible to reach instability at nominally stable conditions given the strong elastic contrast between ice and rock beds (Rice and others, Reference Rice, Lapusta and Ranjith2001). This highlights the range of factors that contribute to seismic instability, and why further experiments and analysis are needed to fully map the conditional dependence of stability.

4. Data processing and machine learning analysis

Each experiment starts with an elastic loading ramp, after which we observe abrupt and audible stress drop events (Supplementary Video 1 and Figs 1b and 2). Most events directly correspond to bulk mechanical stress drops (see Fig. 1b), but to remove AEs associated with other types of sources (such as smaller patches of slip), noisy events, nonevents triggered by noise, and to normalize the waveforms in a way that focuses on the initial wave arrivals (removing secondary wave arrivals), we implemented a data cleaning and normalization approach based on that implemented by Nolte and Pyrak-Nolte (Reference Nolte and Pyrak-Nolte2022).

First, waveforms were trimmed to a total of 1200 data points, including 400 sample points before the trigger point, giving a total window of 120 microseconds. Waveforms were then normalized by the sum of the squared amplitudes of the first 400 sample points after the trigger, multiplied by a cosine taper (effectively weighting the earlier data points using a cosine shape, see Figure S2 in the supplement of Nolte and Pyrak-Nolte [Reference Nolte and Pyrak-Nolte2022]). Zero and large amplitude waveforms were removed, defined as having a sum of the first 400 normalized sample points greater than 15. This threshold was found to give the best catalog of nonnoise events without greatly diminishing the catalog. 325 events were then removed that have high amplitude low frequency noise component. Finally, the waveforms were realigned to the first maximum peak after the trigger, which refined alignment by a few samples or less in most cases. From this catalog of normalized, filtered, and aligned 1200-sample point waveforms, we used a trial-and-error approach to determine how much of the pre- and post-trigger waveforms to use for training the models and found a total length of 150 sample points, with 45 before the trigger, was optimal, giving a 15-microsecond window. This subsample of the waveforms emphasizes the first arrivals of each AE, which are more dependent on source effects, while ignoring the coda (later wave arrivals), which depend more on path effects. Although the original, unprocessed catalog was able to produce as high prediction accuracies (supplementary text S2), not surprisingly since they contain more information, the processed waveforms were clearer to interpret, which is the main point of this study.

After removing noisy waveforms, we end up with 2817 total events, including 1547 waveforms from 6 till experiments and 1270 waveforms from 6 rock experiments, relatively balanced between bed types. With this labeled catalog (Fig. 3), we systematically explored the ability of numerous supervised classification ML algorithms to predict the bed type for each event based on their waveform and spectra. Classification ML models take a training dataset made up of feature vectors, each of which is labeled with a particular ‘ground truth’ class, and the model's parameters are optimized to predict the class of new, previously unseen feature vectors. In our case the features are the normalized amplitudes of each individual sample point in the event's AE waveform or log spectra, and the labels are the origin of the waveform, either a till or rock-bed. Each classification algorithm uses a different mathematical model to map the event to a specific class. Since our analysis focuses on results for the Random Forest Classifier model (Breiman, Reference Breiman2001) we briefly describe the underlying architecture of random forests.

Figure 3. (a) Waveforms plotted in chronological order along y-axis, colored by (normalized) amplitude (red is positive and blue negative). Rock events are plotted on the left and till on the right. (b) Waveforms plotted together for each experiment (labeled on upper left). Sensors exhibit resonance, with waves at the resonant frequency present throughout the recording, even before the arrival, and thus should not affect the prediction. Each waveform (rock in red and till in teal) is plotted with a thin line, so the darker parts show many waveforms aligned on top of each other, and broader lines show less alignment. Since experiments vary significantly by number of events (94–465), that also contributes to the plot of each experiment's appearance. Number of events and bed temperature for each experiment, as well as experiment by experiment training and testing, in order to discount the possibility that experimental differences are being used in the prediction, are explored in supplementary text S3. Although there are subtle visual differences, it is not obvious that the two beds can be deciphered, making it a useful dataset to explore ML-based classification.

A random forest is an ensemble learning model that incorporates the results of multiple decision trees. Each decision tree is a hierarchical, branching structure that splits a dataset into subgroups at each of several decision nodes based on the value of specific features. Random forests average the results of multiple trees, each of which is trained using a random sample of the data (bootstrap aggregating, or bagging) and a random sample of the feature space (feature bagging). This bagging procedure reduces overfitting the model on the training data. Since each decision is made based on a feature of the data, in this case the normalized amplitude at a certain time in the time series or power at a certain frequency in the spectrum of each event, we can inspect the tree to identify which parts of the waveform or which frequencies are most predictive. These ‘feature importances’ therefore allow us to interpret the model from a physical perspective.

In order to train and test the ML algorithms, we choose the most interpretable input features, the normalized waveform amplitudes at each time step or the log₁₀ power at each frequency for the spectra, using the entire 15 microsecond time series. The waveforms and spectra were independently broken into train and test datasets and randomly assigned, containing 70 and 30% of the events, respectively. The trained models select the most important temporal portions of the waveforms or frequencies in the spectra for discriminating between bed labels. We tested a suite of six basic ML algorithms commonly used for classification problems. They are available with the scikit-learn python package (see Open Research section for github page with Jupyter notebook for all processing and analysis steps, as well as supplementary text S2 for more details). The algorithms we tested (and the mean prediction accuracy of each) are Naïve Bayes (~66%), XGBoost (~74%), support vector machines (~77%), random forests (~77%), fully connected neural networks (~75%) and K-nearest Neighbors (~75%). Hyperparameters were tuned for each algorithm and input data type (time or frequency domain) using 5-fold cross validation, and the highest-accuracy model for each algorithm was then used for prediction on the test set. The results of all our tests are summarized in supplementary text S2, but here we focus our analysis on the Random Forest Classifier model (Breiman, Reference Breiman2001) applied to the processed catalog, since it obtained some of our highest prediction accuracies, but, most importantly, the algorithm provides the feature importance needed to interpret how the model obtains its results. By showing the weighting of each waveform sample or frequency in making its prediction (Figs 4a and b), feature importance visually highlights the subtle differences between different event sources. The purpose of this study is to understand how bed differences manifest in the resulting emissions, not to find a black-box algorithm which best differentiates them.

Figure 4. (a) Feature importance (black), showing the weighting of each waveform sample to the model prediction, highlights the importance of the initial, post-trigger, wave arrivals. The superimposed normalized waveforms show till (teal) events are higher amplitude than rock (red) in these first oscillations. (b) Feature importance (black) of each frequency in the model prediction, show till (teal) and rock (red) spectra partially separate from each other above about 1 MHz, with till having more energy at these higher frequencies. It is not clear why the model finds certain frequencies more important for prediction. (c) Distribution of largest repeated mechanical stress-drop amplitude from 23 till and 22 rock experiments at similar conditions show till has overall higher stress-drops, although the two populations overlap significantly. (d) Stress-drops vs recurrence interval for till and rock experiments shows till's greater healing (higher slope) contributes to higher stress-drops, while rock healing varies more, but is generally lower.

5. Bed type classification from acoustic emissions

Using a wide range of classification algorithms, we consistently find prediction accuracy above 50%, mostly between 65 and 80% (supplementary Figure S3), showing it is possible to tell if a population of AEs was emitted by a till or rock-bed. This is not clear by visually examining the waveforms (Fig. 3), so the algorithms successfully extract subtle waveform and spectral features corresponding to the different bed types.

To be able to apply our findings from laboratory AEs to field-scale seismicity, it is vital that we can interpret how the algorithms make their prediction. Although transfer learning methods offer the potential to train with labeled laboratory or modeled datasets and ‘transfer’ the model to more limited field or laboratory data (e.g., Wang and others, Reference Wang, Johnson, Bennett and Johnson2021), clear differences in the spectral content, travel path effects, and scale of field seismic data make this a difficult task. By isolating and interpreting the features the algorithms are using to make their successful predictions, we can identify and quantify predictive features to closely examine for in field data. The feature importance for the Random Forest Classifier model shows that it focuses on the peak and valley of the first full oscillation of the initial wave arrival (Fig. 4a). Plotting all the normalized waveforms (color coded by bed type) together, we can see that the till (teal) waves tend to have higher amplitude in these first peaks. Similarly, log spectra show more energy at higher frequencies for the till in comparison to rock spectra (Fig. 4b). Analyzing the mechanical data from 23 till and 22 rock experiments (including experiments where AEs were not recorded), we find that the stress-drops of stick–slip events on till beds are generally higher (Fig. 4c). The more impulsive arrivals and higher frequency content is consistent with till's higher stress-drops, since seismological stress-drop is calculated by the corner frequency where energy starts to fall off (e.g., Zoet and others, Reference Zoet, Anandakrishnan, Alley, Nyblade and Wiens2012). This, in turn, can be explained by till's higher healing (Figs 2 and 4d), friction (Fig. 2), as well as the rougher till surface (with its larger grain sizes).

It is likely that obtaining much higher prediction accuracies is impossible given the similarity between waveforms. The stress-drop and healing rates of the two populations clearly overlap (Fig. 4c and d); spectra and waveform characteristics do as well. How these similarities affect prediction can be most clearly seen with the log spectra since the visual separation is clearest (Fig. 4b). Figure 5a and b show that misclassified events are in the region between the event types, while Figure 5c shows that the waveform statistical attributes also greatly overlap. Although correctly predicting every event is unrealistic, given a sufficient sample size, our results suggest it could be possible to predict the bed type of a group of events from the same source conditions (see supplementary text S4).

Figure 5. Log spectra of correct and misclassified (a) till and (b) rock events and (c) distributions of statistical measures of waveforms from all experiments from each bed show how much the event populations overlap. The higher variance in the till waveform distributions is due to their more impulsive nature, but there are many rock events with just as high variance.