# Rock slope failure preparation paced by total crack boundary length

Sophie Lagarde, Michael Dietze, Conny Hammer, Martin Zeckra, Anne Voigtländer, Luc Illien, Anne Schöpa, Jacob Hirschberg, Arnaud Burtin, Niels Hovius, et al.

► To cite this version:

Sophie Lagarde, Michael Dietze, Conny Hammer, Martin Zeckra, Anne Voigtländer, et al. Rock slope failure preparation paced by total crack boundary length. Communications Earth & Environment, 2023, 4, 10.1038/s43247-023-00851-0. insu-04155667

HAL Id: insu-04155667

https://insu.hal.science/insu-04155667

Submitted on 7 Jul 2023

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Rock slope failure preparation paced by total crack boundary length

Sophie Lagarde ID 1,2✉, Michael Dietze<fnref n="1" />, Conny Hammer³, Martin Zeckra ID ⁴, Anne Voigtländer ID <fnref n="1" />, Luc Illien<fnref n="1" />,2, Anne Schöpa<fnref n="1" />, Jacob Hirschberg⁵,6, Arnaud Burtin⁷, Niels Hovius<fnref n="1" />,2 & Jens M. Turowski<fnref n="1" />

Gravitational mass wasting prediction requires understanding of the factors controlling failure. Prior to slope failure, cracks in the weakened rock are thought to grow and coalesce, eventually forming a continuous failure plane. Here, we apply a hidden Markov machine learning model to seismic data, revealing the temporal evolution of cracks prior to a major rockslide event in the Swiss Alps. After prolonged linear increase of the crack cumulative number, an S-shaped crack rate pattern occurred in the day before the rockslide. A simple mechanistic model can explain this behaviour, showing that total crack boundary length is the key factor controlling failure plane evolution immediately before mass movement. Our findings imply that cracks should be treated as 2-D, rather than 1-D objects, and that slope failure can be driven predominantly by internal rather than external processes. Our model offers a novel, physically based approach for early warning of slope failures.

<footnote n="1">¹Section 4.6 Geomorphology, GFZ German Research Centre for Geosciences, Potsdam, Germany. ²Institute of Geosciences, University of Potsdam, Potsdam, Germany. ³Institute of Geophysics, University of Hamburg, Hamburg, Germany. ⁴Royal Observatory of Belgium, Uccle, Belgium. ⁵Swiss Federal Institute for Forest, Snow and Landscape Research WSL, Birmensdorf, Switzerland. ⁶Institute of Environmental Engineering, ETH Zurich, Zurich, Switzerland. ⁷Université de Paris, Institut de physique du globe de Paris, CNRS, F-75005 Paris, France. ✉email: lagarde@gfz-potsdam.de</footnote>
Slope failures and landslides are frequent hazard processes in the mountains. They occur suddenly and rapidly and can have high-energy impacts in their transit and deposition areas. As a result, slope failures are responsible for thousands of fatalities each year¹,², as well as billions of dollars of economic damage³,⁴. In addition, they are a dominant process shaping steep topography⁵; they feed the sediment cascade⁶,⁷ and drive biogeochemical cycles⁸,⁹.

Rock slope failure can be conceptualised as a two-phase process¹⁰. In the preparation phase, a volume of rock is progressively detached from the surrounding rock mass by the emergence, growth and coalescence of cracks. This leads to the formation of a failure plane. Once fully detached, an activation phase occurs in which the mass slides, topples or falls downslope. The activation phase can be initiated by an external trigger, such as an earthquake or rainstorm¹¹, but it can also start without a discernible prompt. While the activation phase is well studied¹²⁻¹⁶, considerably less is known about the preparation phase, although it determines the timing, volume and mobilisation mechanism of the failure¹⁷. The inability to observe the state and rate of change of the potential failure planes in natural hillslopes hinders reliable anticipation of slope failure.

Rock cracking emits seismic and acoustic signals that can be detected in the laboratory¹⁸ and in the field¹⁹,²⁰. Thus, a network of geophysical sensors installed at the surface of an unstable rock mass allows indirect yet highly resolved information on processes deep inside the progressively weakening substrate²¹⁻²⁴. However, crack signals are not always intense enough to discern them amidst ambient seismic noise.

Here, we take advantage of the proximity of a series of rock slope failures in the Illgraben catchment, Switzerland, to three instruments in a catchment-wide seismometer network²⁵, with a sufficiently high signal-to-noise ratio (SNR) to allow crack detection. To optimise the seismic detection of cracking in rock slopes and to generate comprehensive data sets necessary to probe the formation of failure planes, we applied a machine learning technique based on hidden Markov models²⁶. We match the observed crack rate evolution prior to a large slope failure with a simple mechanistic model and identify a new first-order control on crack propagation and failure plane evolution prior to mass-wasting.

**A rockslide at the Illgraben, Switzerland.** The Illgraben is a 10 km² steep side valley of the Rhône River²⁷ (Fig. 1). This active debris flow catchment with almost 1900 m local relief²⁸ supplies 5–15% of the total sediment load of the Rhône River upstream of Lake Geneva. The 30–80° steep catchment slopes host frequent rock falls and slides²⁵. The mean annual precipitation is 900 mm yr⁻¹, and the mean annual temperature is 5.9 °C at the basin mean elevation²⁹.

On 2 January 2013, a 10,000 m³ large rockslide happened in the eastern Illgraben flank (860–1560 m in elevation). The preparatory phase and aftermath of this event have been recorded by a network of nine broadband seismometers. Three instruments were located in the vicinity of the event site, one only 180 m from the detachment area (STA6 on Fig. 1). The proximity and positioning of these instruments allowed a systematic recording even of small seismic sources within the unstable slope, providing a rare insight into the rock mass cracking that lead to slope failure.

## Results and discussion

**Crack activity evolution.** Commonly, the onset of seismic events is determined from a sharp rise in the seismic amplitude, for example, by computing the ratio of short-time and long-time average energy and defining a threshold value³⁰. Cracks are

characterised by short (up to a few seconds) seismic signals with mostly 1–20 Hz frequency content that rapidly emerge from the background and fade more slowly¹⁹,²⁰ (Fig. S1a, b). Because many crack signals barely rise above background noise, we applied a machine-learning method based on hidden Markov models (Method), which has been successfully used before to detect events of various types in seismic data²⁶,³¹,³². Our algorithm performance assessment using synthetic data (Method and Supplementary Note 1) indicates that a trend in crack density evolution can reliably be retrieved (Supplementary Note 1 and Supplementary Fig. S2 and S3). We identified 1592 crack events within the 20 days before and 10 days after the January 2nd 2013, main rock slope failure for six seismic stations in the Illgraben network (Fig. 1). In addition, using the same method, we identified slope failure and debris mobilisation during this interval. Slope failures feature an initial rupture followed by a complex fall or slide process that causes rock mass disintegration and downslope movement. The initial rupture distinguishes them from smaller debris mobilisation events, in which material either detaches or becomes remobilized from a prior slope failure deposit.

There was no clear external trigger for the 2013.01.02 rockslide: neither precipitation nor earthquake activity occurred in the five days prior to it. The air temperature crossed the freezing point multiple times within 20 days before the event (Supplementary Fig. S4), implying that snowmelt may have played a role in causing the rockslide but registered at −2 °C at the time of the rockslide. The numbers of detected slope failure and debris mobilisation events showed a sharp increase at the time of the rockslide (Fig. 2a). A second step increase in debris mobilisation activity occurred two days later, soon after the air temperature rose above 5 °C. Notably, the first rise in the number of debris mobilisation events coincided with a change in the trend of the cumulative number of cracks N(t) detected at STA4-6. Whereas N(t) can be approximated by a linear function (R² = 0.989) in the two weeks prior to the rockslide, it displayed an S-shaped development in the last 27 h before failure (green square, Fig. 2a).

**Modelling the crack temporal distribution.** For a slope failure to occur, the complete failure plane area needs to be disjoined (Fig. 3). This is achieved by the growth and coalescence of cracks along an initially largely intact plane. To model this process, we make three simple assumptions:

1. On the failure plane, both the fractured area, $A_f(t)$, and the cumulative number of cracks, $N(t)$, can only increase with time. Given that each crack event creates new fracture area, $N(t)$ scales with $A_f(t)$, and we assume proportionality: $N(t) \propto A_f(t)$. The relationship can be rewritten as $N(t) \propto A(t)$, where $A(t) = A_f(t)/A_{tot}$ is the proportion of fractured area, which increases from zero to one, $A_{tot}$ being the total area of the failure plane.

2. The total crack boundary length $l$ is the sum of the boundary lengths of individual cracks over the whole failure plane. The rate of change of the fractured area scales with the product of $l$ and the crack propagation velocity $v$: $dA_f/dt \sim lv$. We assume that $v$ is constant. Then, $l$ is the only variable controlling the temporal evolution of the fractured area.

3. As the failure plane is initially largely intact ($t = t_0$), the total crack boundary length $l(t_0) \approx 0$ for $A_f(t_0) \approx 0$ (Fig. 3). During the first stages of the failure plane fracturing phase, the number of cracks increases and cracks propagate, leading to an increase in $l(t)$. When a sufficiently large area of the plane is covered by cracks, the cracks begin to join. In this phase, $A_f(t)$ continues to increase, but the behaviour of $l(t)$ is
<figure><img src="image_3.png" /><figcaption>**Fig. 1** Catchment map with the location of 2013.01.02 rockslide and positions of seismic stations (STA). Seismically detected events used for this study happened within the spatial filter area (Method). Inset: Map of Switzerland with the catchment location shown by the white square.</figcaption></figure>

unclear (Fig. 3). Finally, shortly before the failure time $t_f$, when the fractured area approaches the total area, the failure plane is only held by a small number of rock bridges: $l(t \to t_f) \sim 0$ and $A_f(t \to t_f) \sim A_{tot}$ (Fig. 3). As $A_f(t)$ continues to increase, the rock bridge area and consequently $l(t)$ decrease. We use the simplest possible function that is consistent with these boundary conditions, a parabola of the form: $l(A) \sim A_f(1 - A_f)/\tau$, where $\tau$ is a constant.

By combining these equations, we obtained a differential equation (Method) with the solution:

$$N(t) \sim \frac{A_{tot}}{1 + \exp(b - \frac{t}{\tau})} \quad (1)$$

Here, b is the integration constant.
<figure><img src="image_4.png" /><figcaption>**Fig. 2** Cumulative number of geomorphic events as a function of time. **a** Cumulative number of geomorphic events as a function of time from 2012.12.13 to 2013.01.10, retrieved from the seismic records of STA4, STA5 and STA6 using hidden Markov model event detection. The green square highlights the cracks in the 27 h prior to the rockslide shown in (b). The black dashed line corresponds to the timing of the rockslide. **b** Zoom on the cumulative number of cracks displayed in the green square in (a); normalised crack cumulative number per hour for the 27 h prior to the rockslide (green line), and corresponding best fit using equation 1 (black line, τ = 0.1, A<sub>tot</sub> = 0.97, b = 4.3, R<sup>2</sup> = 0.998).</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Fig. 3** Failure plane fracturing as a function of time. As the fractured area $A_f(t)$ (white) grows and the rock bridge area (black) decreases, the total crack boundary length $l(t)$ (green) increases at first and then decreases towards failure.</figcaption></figure>

Equation 1 describes a sigmoidal exponential function, which was fitted to the normalised cumulative number of cracks $N(t)$ for the 27 h period prior to failure, giving $R^2 = 0.998$ (Fig. 2b). The good fit of the model indicates that the total crack boundary

length $l$, an internal state parameter, is the main factor controlling crack evolution in the hours before the event. Other factors, like the crack propagation velocity, which is assumed to be modulated by external controls such as the temperature and humidity$^{33,34}$,
are of secondary importance and can be neglected at first order as they do not improve the model fit (Method, Supplementary Note 2 and Supplementary Fig. S5).

**Two-phase failure evolution from distributed to localised effects.** The proximity of three of our seismometers to the 2013.01.02 rockslide allowed us to distil a unique high-resolution data set on crack occurrence using hidden Markov model event detection. Similar to what has been observed in tunnel construction³⁵ and in small-scale laboratory experiments³⁶, our data suggest the existence of a switch between distributed cracking and localised damage accumulation. In the initial part of our 30-day record, cracks appear stochastically with a constant rate in time (Fig. 2a) and are distributed in space within the entire rock volume (crack initiation phase). This is likely due to external changes in the stress field by chemical and mechanical action, for example, due to water infiltration, temperature changes, or root growth³⁷,³⁸. Each cracking event, in turn, slightly changes the internal stress field of the slope, with stress either decreasing or increasing depending on the location of the considered point with respect to the crack tip³⁹. At some point, crack activity migrates to appear dominantly at the future failure plane. When a threshold in crack density is reached, cracks begin to interact with each other (crack damage phase). At this point, the crack appearance and growth become highly non-linear in time (Fig. 2b) and space, directed by internal constraints, and a failure plane starts to develop³⁵,³⁶.

**Introducing a first-order two-dimensional failure control parameter.** Our study sheds light on the parameters controlling failure plane growth at a slope scale under natural conditions in the hours immediately before an event. The transition from the crack initiation phase to the crack damage phase corresponds to a switch from an external to an internal control on crack expansion. During the crack damage phase, the total crack boundary length becomes the dominant control parameter of failure plane evolution. In previous studies, cracks have usually been treated as one-dimensional objects⁴⁰⁻⁴⁴, and the total crack boundary length was not explicitly considered. Our analysis reveals the importance of the internal state parameter *l* for crack propagation. Our results indicate that a two-dimensional treatment, considering cracks as an areal feature and explicitly including a description of the evolution of the total crack boundary length, is necessary to capture the failure dynamics of rocks at the slope scale.

**Failure events without evident external triggers.** Our analysis suggests that in the hours to days before a major slope failure, the dynamics of cracking are internally controlled. This could explain why it has often proved difficult to identify clear trigger mechanisms for slope failure⁴⁵⁻⁴⁷. Triggers of mass wasting are often condensed into two external factors: seismic shaking and precipitation⁵,⁴⁸. When neither of these triggers is obvious, then attention may turn to a mixture of external parameters, even if none of them shows an obvious connection to the observed slope failure⁴⁹,⁵⁰. These external parameters can, of course, contribute to failure, and we suggest that their role is in modulating the slope's subcritical stress state. Importantly, more care needs to be given to distinguishing between triggering factors responsible for crack damage and the controlling factor contributing to crack growth, keeping in mind that the triggering factor could be an internal parameter rather than an external one. Slope failures without evident external triggers illustrate the relevance of this distinction.

**Pacing failure plane development.** To anticipate mass wasting events and thus construct effective early warning systems, it is imperative to identify the moment in which the slope switches from the preparatory phase to imminent failure. In our model, the parameter sets the time scale of failure plane development, accounting for the steepness of the sigmoidal exponential function (Eq. 1). A reduction in the value of $\tau$ implies a reduction in the time needed to develop the failure plane once the internal mechanism is switched on. For the Illgraben example, we obtained $\tau$ values ranging from 0.1 to 0.15 (Fig. 2b, supplementary Fig. S5), corresponding to 2.7–4 hours. Considering the entire dynamic, this permits a potential maximal warning time of about one day. In landslide-prone settings, early warning systems might combine seismic crack rate monitoring with hidden Markov models for automated crack detection. The onset of a sigmoidal phase in the cumulative crack rate development could be used as an early indicator of an impending slope failure. If the small $\tau$ value retrieved from our Illgraben example is representative, then this sigmoidal phase may last only hours, possibly a day, highlighting the need for real-time monitoring, including crack detection. Further work should focus on the controls on this timescale, for example, by local lithology or climate.

**Method**

**Seismic instrumentation.** From 2012 to 2014, a network of up to ten Nanometrics Trillium Compact 120 s broadband seismometers, sampled by Digos DataCube³ext loggers at 200 Hz, was deployed in and around the catchment of the Illgraben, Vallais, Switzerland, to monitor distributed geomorphic activity²⁵ (Fig. 1). Three broadband seismometers were deployed within a few hundred metres from the location of a major slope failure that mobilised 10,000 m³ of rock. Two stations were located at the base of the slope, to the southwest (STA4) and north (STA5) of the runout zone, and one just 180 m from the edge of the failure scarp (STA6).

**Hidden Markov models workflow.** All computations for the event classification were done in Python 3.8, using the modules Obspy⁵¹ and ASESS (Advanced Seismic Event Spotting System²⁶). ASESS works in three principal steps that have to be repeated for each seismic station: (i) parametrisation of the raw waveform, (ii) model training, and (iii) event detection.

When using hidden Markov models-based classification systems, the detection and classification of signals are done in a single step. Thus, each time step of the incoming continuous data stream is assigned to a specific signal class, which can be either one of the defined event types or noise. As the noise model also captures signals that are not of interest, we call it the background model in the following.

**Parameter value choice.** During the parametrisation step, thirty different features related to spectral, complex trace and polarisation properties of a seismic signal (Table 1 in²⁶) were computed from the waveform in a sliding window fashion. Each window contains 256 samples and is shifted by 1 s over the complete data stream. The use of all 30 available features can lead to redundant information and biased models. To account for that dimensionality problem, we used principal component analysis⁵². Each feature is represented by one eigenvalue and one eigenvector, and the features are sorted by their degree of information. The eigenvalues are comparable to a weighting factor for each dimension. We conducted a principal component analysis on each type of reference event (details below) used to train the model. For each of them, the same two features (normalised envelope and its derivative) had the smallest weight, indicating that they are negligible. Thus, we conducted this analysis on only 28 features. Except for the polarization features, we only used the east component for the computation, as the amplitudes of the reference events (see below) were the highest in this direction. On the main event day, 2 January 2013, seismic energy was maximal within the 1–30 Hz band, typical for mass wasting processes sensed at less than a kilometre distance⁵³. Thus, we applied a bandpass filter to the data on this frequency range. The minimal frequency of 1 Hz discards ocean microseismic effects⁵⁴ and reduces the impact of teleseismic earthquakes⁵⁵.

**Reference events for event classes.** A reference event is an event that is considered representative of a given event class. The definition of reference events is required in steps (i) and (ii) of the hidden Markov models workflow. After a visual inspection of the data, we defined four event classes. Three of them correspond to geomorphological processes: crack, slope failure and debris mobilisation. The crack class is characterised by short (up to a few seconds) seismic signals with mostly 1–20 Hz frequency content that rapidly emerge from the background and fade more slowly¹⁹,²⁰,⁵⁶ (Supplementary Fig. S1a, b). The slope failure class comprises an initial detachment and subsequent downslope transport of the disintegrating
rock mass. It exhibits a short duration of seconds to minutes, with a sharp peak followed by an extended composition of lower amplitude waveforms in the time domain (Supplementary Fig. S1c) and a high energy content in the frequency domain (Supplementary Fig. S1d). The debris mobilisation class, corresponding to the mobilisation or remobilisation of loose slope material without an initial detachment phase, exhibits a short duration of seconds to minutes with emergent waveforms of roughly constant amplitude in the time domain (Supplementary Fig. S1c) and with a lower energy content than for the slope failure class in the frequency domain (Supplementary Fig. S1d). The distinction between these two last classes is justified as they are generated by two different processes. As the Illgraben is located in an active orogen zone with seismic activity, we also defined a tectonic class for local earthquakes, where the first arrival seismic waves are direct crustal phases. Earthquakes are not the subject of our study, but local earthquakes feature a duration and frequency range similar to the geomorphic events described above. Thus, they were included in the classification system in order to prevent misclassification effects of the actual target events.

**Background and event models.** A robust classification system requires a background model that captures all signals that are not of interest. Furthermore, the specific approach implemented in ASESS derives the event model directly from such a background model. Thus, training data for the background class should capture the whole feature space. The geomorphic reference events (crack, slope failure and debris mobilisation) were trained together as they happened within the same hour (2013.01.02, 10:00:00–11:00:00 UTC, Supplementary Fig. S1), with this whole hour used to train the temporary background model. The earthquake reference event was trained separately with the hour of signal around the event (2013.01.05, 07:00:00–08:00:00 UTC) used to train the temporary background model. Finally, the event models were kept for the classification procedure (step iii). Background noise changes over time scales ranging from hours (e.g., temperature variations, anthropogenic activity) days and weeks (e.g., low-pressure areas) to months due to seasonality (e.g., microseisms). If the noise parametrised by the background model differs too much from the noise conditions of the day when the classification is done, it can cause misclassifications. Thus, for every day, a new background model was estimated, and the whole 24 h long record was used for training. The information from the prior model was completely omitted. To avoid bias arising from using the same background model for training and classification, we combined the initial event models with the background model from the preceding day (step iii). In this way, the training and classification data sets differ, but noise variability between both sets is limited. For each day, we obtained a list of the event type, their start time, their duration, and their log-normal likelihood.

**Classification post-processing.** In order to increase the quality of the classification results, we refined them. First, to discard unlikely events, a uniform log-normal likelihood threshold of −70 was applied. Second, based on manually detected events, cracks usually last between 0.5 and 5 s and slope failure and debris mobilisation events last more than 5 s. Thus, we applied these values as event duration thresholds. Third, we took advantage of the network of seismic stations to increase result quality by imposing the events to be detected by the station STA6, located just above the detachment location of the main slope failure zone as well as at least on one at the bottom of the rockslide (STA4 and STA5). All detections not fulfilling these constraints were rejected.

In order to exclude the events that did not happen in the same area (spatial filter, brown box, Fig. 1) as the main event (03:42 UTC, 2013.01.02), the crack events obtained using ASESS were located. The computation was done in Rv4.2.0⁵⁷, using the “eseis” package v. 0.6.0⁵⁸. To locate the events, at least three seismic stations are needed. Regardless of the events detected by the machine learning algorithm, any event generates seismic waves travelling in every direction. Because the stations STA6, STA5 and STA4 were located in the same area, it is reasonable to expect that one event detected by at least two of them is also recorded on the third, even though the machine learning algorithm did not recognise it. Thus, in order to maximise the number of crack events with low SNR, we included events detected on at least two of the three stations closest to the rockslide, i.e. STA4, STA5 and STA6. Afterwards, all three stations were used to locate the events. The seismic data were deconvolved, detrended and filtered between 1 and 30 Hz. The envelope of the signal was computed. We used a signal migration approach²⁵. In this method, empirical envelope time delays among station pairs are estimated by cross-correlation analysis. Distance maps corresponding to grids for each seismic station containing the distance of a pixel to the station are generated. We used a DEM with 30 m resolution⁵⁷ to correct for topographic effects. The distance maps were converted to time delay tables using an apparent seismic wave velocity (here 500 m s⁻¹, from²⁵). The empirical time delays are then compared to the synthetic times, giving the theoretical location of the seismic source as a grid with values ranging between 0 and 1, denoting the probability density of the location estimate. Then, a density map estimating the point distribution in space was done using the function “ppp”, from the “spatstat” package v. 2.3-0⁵⁸, with a value of 10 for the standard deviation.

The events that occurred outside the area of the main event (03:42 UTC, 2013.01.02) were excluded from the subsequent analysis (spatial filter, brown box, Fig. 1). Assuming that cracks not directly located on the failure plane also

contribute to stress release and considering the uncertainty of the location algorithm, the spatial filter was designed somewhat wider than the rockslide area (Supplementary Fig. S7).

## Hidden Markov model performance assessment

**Synthetic testing.** AESS has the advantage of providing an automatic classification. The use of the feature domain associated with hidden Markov models would, in principle, allow for retrieval of signals with low SNR that a user would have probably missed. However, AESS can also give false positive results. As crack events have low SNR, we conducted a synthetic data study to estimate the confidence level of the classification.

To explore the impact of SNR on crack detection, we created four synthetic data sets. All computations were done in R. The basic idea was to use an empiric crack signal, convolve noise to it in the frequency domain, and investigate the detection reliability of the model. We conducted all of the tests computing the same features that were used to process the Illgraben data, excluding the polarisation features, i.e., only 22 features were kept in the parameterisation. Conducting the tests using the polarization features would have meant generating noise on the three different channels. However, even if noise is random by definition, we cannot safely exclude between-channel correlation effects, which would be hard to propagate into synthetic data sets.

Two types of noise were applied: white noise and pink noise. White noise was tested as it is the simplest kind of noise. Some empirically measured noise spectra of the Illgraben, i.e., time series for which no event was detected during the classification, revealed that the noise encountered there is actually pink noise. Thus, we fitted the empirically measured noise spectra of the Illgraben with a power scale function $f(x) = 10(a^*\log_{10}(x)+b)$ using a nonlinear least square regression function. We then added white noise to this initial function in order to obtain pink noise.

For both types of noise, we created two different test data sets. For the first test data set, the crack event used for the model training and the classification is the same; only the noise signal differs. For the second test data set, the crack event used for the model training and the noise signal are different from the ones used in the classification. For each test data set, the amplitude of the noise was varied in order to obtain SNR ranging from 1 to above 9.

For each test, the SNR was computed. Considering an event of length *d*, with noise<sub>pre</sub> the noise right before the event, having a length *d*/2, and noise<sub>post</sub> the noise right after the event, having a length *d*; we define the SNR as the ratio of the maximum of the signal envelope and the mean of the envelope of noise<sub>pre</sub> and noise<sub>post</sub> aggregated together. To allow comparison, the SNR of the cracks detected on the Illgraben dataset using AESS were also computed following this definition (Supplementary Fig. S3, a), as well as the temporal crack distribution evolution per SNR (Supplementary Fig. S3, b).

We define a crack as falsely detected (false positive) if its starting time is detected by AESS after the end of the crack or if it was over before the start of the crack. We define a crack as correctly detected (true positive) if its detected time span overlaps partly with the empiric crack. For each SNR, we computed the percentage of false positives as well as true positives (Supplementary Fig. S2).

**Confusion matrix.** We computed a confusion matrix to give additional information about the confusion between the different classes. Prior to the application of the hidden Markov models, no catalogue of the events recorded during the studied period (2012.12.13 to 2013.01.10) existed. Establishing such a catalogue and performing a visual examination of the seismic data is time-consuming and out of the scope of this study. Thus, in order to compute the confusion matrix, a visual examination of the seismic trace was performed in a short time period (7 h). From 2013.01.01, 20:00:00 UTC to 2013.01.02, 3:00:00 UTC, 13% of the slope failure events (23 events), 18% of the debris mobilisation events (93 events) and 1.7% of the crack events (29 events) were recorded by the hidden Markov models: this time period offers a good compromise between the number of events recorded and the amount of seismic data to visually assess.

In order to establish the visual catalogue, the seismic trace of the east component of STA6 was visually inspected in the time domain after the application of a bandpass filter between 1 and 30 Hz (same component and bandpass filter applied to the seismic signal prior to the application of the machine learning technique). The visual catalogue was then compared to the list of events obtained after the two first steps of the post-processing (the final list prior to the removal of the events not located in the area of interest), and the confusion matrix was constructed using this comparison.

**Fitting of the physical model to the data.** The physical model was built according to the assumptions given in the main text, paragraph “Modelling the crack temporal distribution”. For the convenience of the reader, we detail here the steps to obtain the final equations for the cumulative crack number *N*(*t*) as well as the non-cumulative temporal crack number *dN*/*dt*. A differential equation is obtained by combining assumptions 2 and 3 from the main text:

$$ \frac{dA_f}{dt} \sim \frac{A_f(1 - A_f)}{\tau} $$
This equation can be integrated to yield:

$$A_f(t) \sim \frac{1}{1 + \exp(b - \frac{t}{\tau})}$$

Taking into account assumption 1, the cumulative crack number $N(t)$ is obtained (Eq. 1, main text). The non-cumulative temporal crack evolution $dN/dt$ corresponds to the derivative of $N(t)$ at time $t$. Its expression is:

$$\frac{dN}{dt} \sim \frac{A_{\text{tot}}}{\tau} \frac{\exp(b - \frac{t}{\tau})}{(1 + \exp(b - \frac{t}{\tau}))^2} \quad (\text{M1})$$

All of the fits between the model and the data were done using the function “curve_fit” from the package “scipy.optimise” in python (Fig. 2 and Supplementary Fig. S5).

**Testing the physical model with a non-constant velocity.** In order to test the impact of the crack propagation velocity formulation on the model results, we replace the assumption that $v$ is a constant, as made in assumption 2, with the expression $v \sim v_0e^{et/t_0}$ (Supplementary Fig. S5). The new expressions obtained for $N(t)$ and $dN/dt$ then become:

$$N(t) \sim \frac{A_{\text{tot}}}{1 + \exp(b - t_0 v_0 e^{et/t_0} / \tau)}$$

$$\frac{dN}{dt} \sim \frac{A_{\text{tot}} v_0 \exp(t_0 v_0 e^{et/t_0} / \tau + b + t/t_0)}{\tau(\exp(t_0 v_0 e^{et/t_0} / \tau) + e^b)^2}$$

To avoid having too many free parameters, we fixed $v_0$ and $t_0$ to 1, giving:

$$N(t) \sim \frac{A_{\text{tot}}}{1 + \exp(b - e^t / \tau)} \quad (\text{M2})$$

$$\frac{dN}{dt} \sim \frac{A_{\text{tot}} \exp(e^t / \tau + b + t)}{\tau(\exp(e^t / \tau) + e^b)^2} \quad (\text{M3})$$

## Data availability

The raw seismic data is available at https://doi.org/10.14470/4W615776⁵⁹.

## Code availability

The method and supplementary information contain the information necessary to reproduce the results of this study. All computations for the event classification were done in Python 3.8, using the modules Obspy⁵¹ and ASESΣ (Advanced Seismic Event Spotting System²⁶). The event location in the post-processing stage was done in R⁶⁰ with the “eseis” package v. 0.6.0¹⁰. The fits between the model and the data were done using the function “curve_fit” from the package “scipy.optimise” in python. All these packages and functions are publicly available in the references.

Received: 19 September 2022; Accepted: 16 May 2023;
Published online: 05 June 2023

## References

1. Froude, M. J. & Petley, D. N. Global fatal landslide occurrence from 2004 to 2016. *Nat. Hazards Earth Syst. Sci.* **18**, 2161–2181 (2018).

2. Petley, D. Global patterns of loss of life from landslides. *Geology* **40**, 927–930 (2012).

3. Dilley, M. *Natural Disaster Hotspots: A Global Risk Analysis* (Vol. 5). (World Bank Publications, 2005).

4. Kjekstad, O. & Highland, L. Economic and social impacts of landslides. in *Landslides–Disaster Risk Reduction* (pp. 573–587). (Springer, Berlin, Heidelberg, 2009).

5. Densmore, A. L. & Hovius, N. Topographic fingerprints of bedrock landslides. *Geology* **28**, 371–374 (2000).

6. Hovius, N., Stark, C. P. & Allen, P. A. Sediment flux from a mountain belt derived by landslide mapping. *Geology* **25**, 231–234 (1997).

7. Korup, O., Densmore, A. L. & Schlunegger, F. The role of landslides in mountain range evolution. *Geomorphology* **120**, 77–90 (2010).

8. Battin, T. J. et al. The boundless carbon cycle. *Nat. Geosci.* **2**, 598–600 (2009).

9. Ruttenberg, K. C. The global phosphorus cycle. *Treatise Geochem.* **8**, 682 (2003).

10. Crozier, M. J. & Glade, T. Landslide hazard and risk: issues, concepts and approach. *Landslide Hazard and Risk.* pp. 1–40 (2005).

11. Wieczorek, G. F. Landslides: investigation and mitigation. Chapter 4— landslide triggering mechanisms. *Transportation Research Board Special Report* (1996).

12. Brodsky, E. E., Gordeev, E. & Kanamori, H. Landslide basal friction as measured by seismic waves. *Geophys. Res. Lett.* **30**, 2236 (2003).

13. Favreau, P., Mangeney, A., Lucas, A., Crosta, G. & Bouchut, F. Numerical modeling of earthquakes. *Geophys. Res. Lett.* **37**, 305 (2010).

14. McDougall, S. 2014 Canadian Geotechnical Colloquium: landslide runout analysis—current practice and challenges. *Can. Geotech. J.* **54**, 605–620 (2017).

15. Moretti, L., et al. Numerical modeling of the Mount Steller landslide flow history and of the generated long period seismic waves. *Geophys. Res. Lett.* **39**, 402 (2012).

16. Schöpa, A. et al. Dynamics of the Askja caldera July 2014 landslide, Iceland, from seismic signal analysis: precursor, motion and aftermath. *Earth Surf. Dyn.* **6**, 467–485 (2018).

17. Chae, B. G., Park, H. J., Catani, F., Simoni, A. & Berti, M. Landslide prediction, monitoring and early warning: a concise review of state-of-the-art. *Geosci. J.* **21**, 1033–1070 (2017).

18. Brantut, N., Baud, P., Heap, M. J. & Meredith, P. G. Micromechanics of brittle creep in rocks. *J. Geophys. Res.* **117**, 412 (2012).

19. Dietze, M., Krautblatter, M., Illien, L. & Hovius, N. Seismic constraints on rock damaging related to a failing mountain peak: the Hochvogel, Allgäu. *Earth Surf. Process. Landf.* **46**, 417–429 (2021).

20. Weber, S., Faillettaz, J., Meyer, M., Beutel, J. & Vieli, A. Acoustic and microseismic characterization in steep bedrock permafrost on Matterhorn (CH). *J. Geophys. Res.* **123**, 1363–1385 (2018).

21. Fiolleau, S., et al. Seismic characterization of a clay-block rupture in Harmalière landslide, French Western Alps. *Geophys. J. Int.* https://doi.org/10.1093/gji/ggaa050 (2020).

22. Amitrano, D., Grasso, J. R. & Senfaute, G. Seismic precursory patterns before a cliff collapse and critical point phenomena. *Geophys. Res. Lett.* **32**, L08314 (2005).

23. Senfaute, G., Duperret, A. & Lawrence, J. Micro-seismic precursory cracks prior to rockfall on coastal chalk cliffs: a case study at Mesnil-Val, Normandie, France. *Nat. Hazards Earth Syst. Sci.* **9**, 1625–1641 (2009).

24. Levy, C., Jongmans, D. & Baillet, L. Analysis of seismic signals recorded on a prone-to-fall rock column (Vercors massif, French Alps). *Geophys. J. Int.* **186**, 296–310 (2011).

25. Burtin, A., Hovius, N., McArdell, B. W., Turowski, J. M. & Vergne, J. Seismic constraints on dynamic links between geomorphic processes and routing of sediment in a steep mountain catchment. *Earth Surf. Dyn.* **2**, 21–33 (2014).

26. Hammer, C., Beyreuther, M. & Ohrnberger, M. A seismic-event spotting system for volcano fast-response systems. *Bull. Seismol. Soc. Am.* **102**, 948–960 (2012).

27. Schlunegger, F. & Norton, K. P. Water versus ice: the competing roles of modern climate and Pleistocene glacial erosion in the Central Alps of Switzerland. *Tectonophysics* **602**, 370–381 (2013).

28. Badoux, A., Graf, C., Rhyner, J., Kuntner, R. & McArdell, B. W. A debris-flow alarm system for the Alpine Illgraben catchment: design and performance. *Nat. Hazards* **49**, 517–539 (2009).

29. Hirschberg, J., Badoux, A., McArdell, B. W., Leonarduzzi, E. & Molnar, P. Evaluating methods for debris-flow prediction based on rainfall in an Alpine catchment. *Nat. Hazards Earth Syst. Sci.* **21**, 2773–2789 (2021).

30. Allen, R. Automatic phase pickers: their present use and future prospects. *Bull. Seismol. Soc. Am.* **72**, S252–S242 (1982).

31. Dammeier, F., Moore, J., Hammer, C., Hasslinger, F. & Loew, S. Automatic detection of alpine rockslides in continuous seismic data using Hidden Markov Models. *Geophys. Res.* **121**, 351–371, https://doi.org/10.1002/2015JF003647 (2016).

32. Niemz, P. et al. Full-waveform-based characterization of acoustic emission activity in a mine-scale experiment: a comparison of conventional and advanced hydraulic fracturing schemes. *Geophys. J. Int.* **222**, 189–206 (2020).

33. Nara, Y., Hiroyoshi, N., Yoneda, T. & Kaneko, K. Effects of relative humidity and temperature on subcritical crack growth in igneous rock. *Int. J. Rock. Mech. Min. Sci.* **47**, 640–646 (2010).

34. Waza, T., Kurita, K. & Mizutani, H. The effect of water on the subcritical crack growth in silicate rocks. *Tectonophysics* **67**, 25–34 (1980).

35. Diederichs, M. S. The 2003 Canadian Geotechnical Colloquium: mechanistic interpretation and practical application of damage and spalling prediction criteria for deep tunnelling. *Can. Geotech. J.* **44**, 1082–1116 (2007).

36. Innocente, J. C., Paraskevopoulou, C. & Diederichs, M. S. Estimating the long-term strength and time-to-failure of brittle rocks from laboratory testing. *Int. J. Rock. Mech. Min. Sci.* **147**, 104900 (2021).

37. Voigtländer, A., Leith, K. & Krautblatter, M. Subcritical crack growth and progressive failure in Carrara marble under wet and dry conditions. *J. Geophys. Res.* **123**, 3780–3798 (2018).

38. Wang, H. F., Bonner, B. P., Carlson, S. R., Kowallis, B. J. & Heard, H. C. Thermal stress cracking in granite. *J. Geophys. Res.* **94**, 1745–1758 (1989).

39. Lin, X., Chu, R. & Zeng, X. Rupture processes and Coulomb stress changes of the 2017 Mw 6.5 Jiuzhaigou and 2013 Mw 6.6 Lushan earthquakes. *Earth Planets Space* **71**, 1–15 (2019).
40. Amitrano, D. & Helmstetter, A. Brittle creep, damage, and time to failure in rocks. *J. Geophys. Res.* **111**, 201 (2006).

41. Baud, P., Reuschlé, T. & Charlez, P. An improved wing crack model for the deformation and failure of rock in compression. In *International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts*. Vol. 33, pp. 539–542. (Pergamon, 1996).

42. Brantut, N., Heap, M. J., Meredith, P. G. & Baud, P. Time-dependent cracking and brittle creep in crustal rocks: a review. *J. Struct. Geol.* **52**, 17–43 (2013).

43. Kemeny, J. M. A model for non-linear rock deformation under compression due to sub-critical crack growth. In *International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts*. Vol. 28, pp. 459–467. (Pergamon, 1991).

44. Main, I. G. A damage mechanics model for power-law creep and earthquake aftershock and foreshock sequences. *Geophys. J. Int.* **142**, 151–161 (2000).

45. Dietze, M. et al. Impact of nested moisture cycles on coastal chalk cliff failure revealed by multiseasonal seismic and topographic surveys. *J. Geophys. Res.* **125**, e2019JF005487 (2020).

46. Helmstetter, A. & Garambois, S. Seismic monitoring of Séchilienne rockslide (French Alps): Analysis of seismic signals and their correlation with rainfalls. *J. Geophys. Res.* **115**, 16 (2010).

47. Stock, G. M., et al. *Historical rock falls in Yosemite National Park, California, 1857–2011*. p. 17. (US Department of the Interior, US Geological Survey, 2013).

48. Meunier, P., Hovius, N. & Haines, J. A. Topographic site effects and the location of earthquake induced landslides. *Earth Planet. Sci. Lett.* **275**, 221–232 (2008).

49. Dietze, M., Turowski, J. M., Cook, K. L. & Hovius, N. Spatiotemporal patterns, triggers and anatomies of seismically detected rockfalls. *Earth Surf. Dyn.* **5**, 757–779 (2017).

50. Guzzetti, F., Peruccacci, S., Rossi, M. & Stark, C. P. Rainfall thresholds for the initiation of landslides in central and southern Europe. *Meteorol. Atmos. Phys.* **98**, 239–267 (2007).

51. Beyreuther, M. et al. ObsPy: A Python toolbox for seismology. *Seismol. Res. Lett.* **81**, 530–533 (2010).

52. Song, F., Guo, Z. & Mei, D. Feature selection using principal component analysis. In *2010 International Conference on System Science, Engineering Design and Manufacturing Informatization*. Vol. 1, pp. 27–30. (IEEE, 2010).

53. Burtin, A., Hovius, N. & Turowski, J. M. Seismic monitoring of torrential and fluvial processes. *Earth Surf. Dyn.* **4**, 285–307 (2016).

54. Longuet-Higgins, M. S. A theory of the origin of microseisms. *Philos. Trans. R. Soc. Lond. Ser. A Math. Phys. Sci.* **243**, 1–35 (1950).

55. Boatwright, J. & Choy, G. L. Teleseismic estimates of the energy radiated by shallow earthquakes. *J. Geophys. Res.* **91**, 2095–2112 (1986).

56. Provost, F. et al. Towards a standard typology of endogenous landslide seismic sources. *Earth Surf. Dyn.* **6**, 1059–1088 (2018).

57. JAXA—Japan Aerospace Exploration Agency. ALOS World 3D 30 meter DEM. V3.2, Jan 2021. Distributed by OpenTopography. https://doi.org/10.5069/G94M92HB. Accessed 14.01.2022 (2021).

58. Baddeley, A., Rubak, E. & Turner, R. *Spatial Point Patterns: Methodology and Applications With R*. (CRC Press, 2015).

59. Burtin, A., et al. Illgraben. GFZ Data Services. *Other/Seismic Network*. https://doi.org/10.14470/4W615776 (2023).

60. R Core Team. R: A language and environment for statistical computing [Computer software manual]. Vienna, Austria. Retrieved from https://www.R-project.org/ (2020).

61. Dietze, M. The R package “eseis”-a software toolbox for environmental seismology. *Earth Surf. Dyn.* **6**, 669–686 (2018).

## Acknowledgements

This work was funded by GFZ internal funding.

## Author contributions

The study is based on M.Z.'s master's research under the supervision of C.H. and N.H., where hidden Markov models were applied to the Illgraben data set. Their results were reproduced, refined and extended by S.L. and J.M.T., who also developed the physical interpretation and model. S.L. and M.D. conducted the synthetic data analysis. J.H. provided and analysed weather data. A.B. was the P.I. of the seismic data collection project at the Illgraben. A.V. made Fig. 3. S.L. and J.M.T. wrote the paper with input from all authors, including A.S. and L.I.

## Funding

Open Access funding enabled and organized by Projekt DEAL.

## Competing interests

The authors declare no competing interests.

## Additional information

**Supplementary information** The online version contains supplementary material available at https://doi.org/10.1038/s43247-023-00851-0.

**Correspondence** and requests for materials should be addressed to Sophie Lagarde.

**Peer review information** *Communications Earth & Environment* thanks Danica Roth, Clément Hibert and the other anonymous reviewer(s) for their contribution to the peer review of this work. Primary Handling Editor: Joe Aslin. Peer reviewer reports are available.

**Reprints and permission information** is available at http://www.nature.com/reprints

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

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

**Open Access** This article is licensed under a Creative Commons

Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2023