Introduction

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 year1,2, as well as billions of dollars of economic damage3,4. In addition, they are a dominant process shaping steep topography5; they feed the sediment cascade6,7 and drive biogeochemical cycles8,9.

Rock slope failure can be conceptualised as a two-phase process10. 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 rainstorm11, but it can also start without a discernible prompt. While the activation phase is well studied12,13,14,15,16, considerably less is known about the preparation phase, although it determines the timing, volume and mobilisation mechanism of the failure17. 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 laboratory18 and in the field19,20. 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 substrate21,22,23,24. 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 network25, 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 models26. 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 km2 steep side valley of the Rhône River27 (Fig. 1). This active debris flow catchment with almost 1900 m local relief28 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 slides25. The mean annual precipitation is 900 mm yr−1, and the mean annual temperature is 5.9 °C at the basin mean elevation29.

Fig. 1: Catchment map with the location of 2013.01.02 rockslide and positions of seismic stations (STA).
figure 1

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.

Full size image

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 value30. 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 slowly19,20 (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 data26,31,32. 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\left(t\right)\) detected at STA4-6. Whereas \(N\left(t\right)\) can be approximated by a linear function (\({R}^{2}=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).

Fig. 2: Cumulative number of geomorphic events as a function of time.
figure 2

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, \(\tau =0.1\), \({A}_{{tot}}=0.97\), \(b=4.3\), \({R}^{2}=\) 0.998).

Full size image

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. 1.

    On the failure plane, both the fractured area, \({A}_{f}(t)\), and the cumulative number of cracks, \(N\left(t\right)\), 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}\left({t}\right)\propto {A}_{f}\left({t}\right)\). The relationship can be rewritten as \({N}\left({t}\right)\propto A\left({t}\right)\), where \(A(t)={A}_{f}(t){/A}_{{{{{\rm{tot}}}}}}\) is the proportion of fractured area, which increases from zero to one, \({A}_{{{{{\rm{tot}}}}}}\) being the total area of the failure plane.

  2. 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. 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 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}_{{{{{\rm{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({{{{{\rm{A}}}}}}) \sim {A}_{f}(1-{A}_{f})/\tau\), where \(\tau\) is a constant.

Fig. 3: Failure plane fracturing as a function of time.
figure 3

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.

Full size image

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

$$N\left(t\right) \sim \frac{{A}_{{tot}}}{1+{{{{{\rm{exp }}}}}}\left(b-\frac{t}{\tau }\right)}$$
(1)

Here, b is the integration constant.

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 R2 = 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 humidity33,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 construction35 and in small-scale laboratory experiments36, 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 growth37,38. 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 tip39. 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 develop35,36.

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 objects40,41,42,43,44, 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 failure45,46,47. Triggers of mass wasting are often condensed into two external factors: seismic shaking and precipitation5,48. 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 failure49,50. 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 activity25 (Fig. 1). Three broadband seismometers were deployed within a few hundred metres from the location of a major slope failure that mobilised 10,000 m3 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 Obspy51 and ASESS (Advanced Seismic Event Spotting System26). 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 in26) 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 analysis52. 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 distance53. Thus, we applied a bandpass filter to the data on this frequency range. The minimal frequency of 1 Hz discards ocean microseismic effects54 and reduces the impact of teleseismic earthquakes55.

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 slowly19,20,56 (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.057, using the “eseis” package v. 0.6.058. 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 approach25. 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 resolution57 to correct for topographic effects. The distance maps were converted to time delay tables using an apparent seismic wave velocity (here 500 m s−1, from25). 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-058, 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

ASESS has the advantage of providing an automatic classification. The use of the feature domain associated with hidden Markov models would, in principle, allow for to retrieval of signals with low SNR that a user would have probably missed. However, ASESS 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*log10(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 noisepre the noise right before the event, having a length d/2, and noisepost the noise right after the event, having a length d/2; we define the SNR as the ratio of the maximum of the signal envelope and the mean of the envelope of noisepre and noisepost aggregated together. To allow comparison, the SNR of the cracks detected on the Illgraben dataset using ASESS 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 ASESS 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{d{A}_{f}}{{dt}} \sim \frac{{A}_{f}\left(1-{A}_{f}\right)}{\tau }$$

This equation can be integrated to yield:

$${A}_{f}\left(t\right) \sim \frac{1}{1+\exp \left(b-\frac{t}{\tau }\right)}$$

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}_{{{{{\rm{tot}}}}}}}{\tau }\frac{\exp \left(b-\frac{t}{\tau }\right)}{{\left(1+\exp \left(b-\frac{t}{\tau }\right)\right)}^{2}}$$
(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}_{0}{e}^{t/{t}_{0}}\) (Supplementary Fig. S5). The new expressions obtained for \(N\left(t\right)\) and \({dN}/{dt}\) then become:

$$N\left(t\right) \sim \frac{{A}_{{tot}}}{1+\exp \left(b-{t}_{0}{v}_{0}{e}^{t/{t}_{0}}/\tau \right)}$$
$$\frac{dN}{dt} \sim \frac{{A}_{{{{{\rm{tot}}}}}}{v}_{0}\exp ({t}_{0}{v}_{0}{e}^{t/{t}_{0}}/\tau +b+t/{t}_{0})}{\tau {(\exp ({t}_{0}{v}_{0}{e}^{t/{t}_{0}}/\tau) +{e}^{b})}^{2}}$$

To avoid having too many free parameters, we fixed \({v}_{0}\) and \({t}_{0}\) to 1, giving:

$$N\left(t\right) \sim \frac{{A}_{{tot}}}{1+\exp \left(b-{e}^{t}/\tau \right)}$$
(M2)
$$\frac{{dN}}{{dt}} \sim \frac{{A}_{{{{{\rm{tot}}}}}}\exp \left({e}^{t}/\tau +b+t\right)}{{\tau \left(\exp \left({e}^{t}/\tau \right)+{e}^{b}\right)}^{2}}$$
(M3)