Théotime De La Selle
email: theotime.de-la-selle@insa-lyon.fr

Stéphanie Deschanel
email: stephanie.deschanel@insa-lyon.fr

Jérome Weiss
email: jerome.weiss@univ-grenoble-alpes.fr

Acoustic emission classification by cross-correlation

Keywords: Acoustic emission, fatigue crack growth, multiplets, cross-correlation, waveform clustering

Acoustic emission signals specific to the incremental advance of fatigue cracks can be detected, cycle after cycle, during crack propagation tests. They were termed acoustic multiplets in reference to repeating earthquakes. Detecting such multiplets would give information about the fracturing process and provide early warnings when monitoring industrial parts in service. We developed a clustering method to extract multiplets signals from all signals and noise of fatigue tests. This method is based on the computation of the cross-correlation function between waveforms. On the distribution of the cross-correlation matrix, signals belonging to multiplets are clearly identified and permit to automatically clustered. This paper developps the steps of the method, presents results obtained on a dataset of a fatigue test and shows its noise robustness.

Introduction

Despite many efforts over the last hundred years to understand its causes and detect its possible precursors, cyclic fatigue of materials remains one of the major causes of failure of industrial parts [START_REF] Suresh | Fatigue of materials[END_REF]. The problem is "How to alert in real time when a crack is likely to lead to the ruin of the part ?". Acoustic emission (AE) is a key method of non-destructive testing, as it enables the recording of transient elastic waves occurring during damage of materials and structures. Numerous works carried out with this method have shown that there is a correlation between the global acoustic activity and damage, including during fatigue [START_REF] Harris | Continuous monitoring of fatigue-crack growth by acoustic emission techniques[END_REF]. However, this correlation only becomes significant near the final failure and thus far too late [START_REF] Berkovits | Study of fatigue-crack characteristics by acoustic emission[END_REF]. Indeed, the monitoring of slow crack growth from a global measurement is difficult and very sensitive to the signal to noise ratio (SNR). Moreover, the non-specific nature of these AE measurements makes the identification of sources difficult [START_REF] Deschanel | Contrôle de la fatigue des matériaux par émission acoustique[END_REF].

A new approach, recently proposed [START_REF] Deschanel | Contrôle de la fatigue des matériaux par émission acoustique[END_REF][START_REF] Deschanel | Acoustic emission multiplets as early warnings of fatigue failure in metallic materials[END_REF], allows the detection of acoustic emissions specific of fatigue crack growth, called "acoustic multiplets" in reference to the analogous phenomenon in seismology [START_REF] Lengliné | Inferring the Coseismic and Postseismic Stress Changes caused by the 2004, M=6 Parkfield Earthquake from Variations of Recurrence Times of Microearthquakes[END_REF]. These so-called acoustic multiplets are characterized by highly correlated waveforms (figure 1 a), signature of a unique source. They are repeatedly triggered over many successive loading cycles at almost the same stress level (figure 1 c, d ande) and originate from a single location (figure 1 a). They are considered to mark the slow, incremental propagation of a crack at each cycle (see striations on figure 1 b), or the rubbing along its faces (fretting). Being specific to incremental fatigue cracking, they can therefore potentially be used as early warnings of global failure of industrial parts [START_REF] Deschanel | Acoustic emission multiplets as early warnings of fatigue failure in metallic materials[END_REF]. A key point in this approach is the choice of an appropriate and robust measure of waveforms similarity. Due to the classical procedure to detect acoustic bursts from threshold crossing and the noise that affect the signal, nearly identical signals can be shifted by a little time delay and are then non-aligned. It means that the method of measure needs to be able to quantify the similarity between 2 waveforms regardless of the delay.

The cross-correlation function is often chosen for this kind of problematic [START_REF] Proakis | Digital Signal Processing[END_REF][START_REF] Kurz | Similarity matrices as a new feature for acoustic emission analysis of concrete[END_REF][START_REF] Grabec | Application of correlation techniques for localization of acoustic emission sources[END_REF]. In the field of acoustic emission, the cross-correlation function has been sometimes used to improve the localisation of acoustic sources [START_REF] Grabec | Application of correlation techniques for localization of acoustic emission sources[END_REF]. A cross-correlation matrix over all the AE signals signals of a channel recorded during a mechanical test on concrete has been used to investigate very quickly changes in the physical properties of the ray path of events from one cluster [START_REF] Kurz | Similarity matrices as a new feature for acoustic emission analysis of concrete[END_REF]. Another application of the cross-correlation matrix has been proposed to filter noise in fatigue test [START_REF] Emamian | Robust Clustering of Acoustic Emission Signals Using Neural Networks and Signal Subspace Projections[END_REF]. This article offers a robust way based on cross-correlation matrix distribution to extract highly correlated populations of signals, in order to identify multiplets in fatigue crack growth tests.

Method to extract highly correlated populations

Cross-correlation function

Cross-correlation consists in the displaced dot product between two signals. It is often used to quantify the degree of similarity or interdependence between two signals [START_REF] Proakis | Digital Signal Processing[END_REF]. In the case of AE, since all measurements were recorded using digital acquisition systems, signals have been discretized, so that the cross-correlation between two signals v and w with the same N samples length is expressed by equation 1.

(1)

When the discrete time series v and w match, the value of corr[v,w] is maximized. We called the maximizing point t m : it corresponds to the time delay between the two waveforms. This is explained when peaks (positive areas) are aligned, making a large contribution to the summation. In the case of AE waveforms, the two parameters n 0 and N define a cross-correlation window. Because acoustic waveforms are transient waves and contain a pre-trigger, an impulsive part and a coda, the choice of this window has to be set precisely for an adequate determination of the maximizing point.

Cross-correlation matrix

By calculating all cross-correlation functions for each pair of waveforms for a single channel, and searching for the maximum of these functions, it is possible to build a crosscorrelation matrix. Hence, the coefficients contain a value in [0, 1] corresponding to the normalized measure of similarity by cross-correlation (see equation 2).

(2

)
An example is presented in fig. 2 showing a cross-correlation matrix of 198 signals compared to each others. By construction, it is a symmetrical matrix with a diagonal full of 1. Then, only half of the matrix is computed.

Fig. 2 Example of cross-correlation matrix calculated from a fatigue crack growth test on steel

Since multiplets are groups of highly similar waveforms emitted close in time, multiplets appear as yellow squares on fig. 2. In the rest of the article, for the sake of simplicity, the term 'crosscorrelation' will refer to the value taken by a coefficient of the cross-correlation matrix as described in equation 2.

Distribution of cross-correlation

An operator can visually identify groups of waveforms highly similar in the crosscorrelation matrix, but an automatized way is possible by calculating the distribution of the matrix. In a fatigue test characterized by acoustic multiplets, two peaks are expected in this distribution (see fig. 3) : one narrow peak close to 1 containing cross-correlations between waveforms belonging to multiplets, and a wider one corresponding to uncorrelated source signals. The remaining correlation of these uncorrelated source signals results from the resonant nature of the AE sensors. Indeed, due to the acquisition chain, and particularly to the frequency domain of the sensor, the cross-correlation between two totally different source signals is larger than zero. Nevertheless, a distribution of cross-correlations computed on a set of data from an experiment without multiplets would not exhibit the peak of high cross-correlations.

Fig. 3 Schematic cross-correlation distribution for a fatigue test emitting multiplets

From this distribution, it is easy to extract automatically multiplets by setting a threshold on the cross-correlation , e.g. 0.8 for the schematic figure above. All waveforms implicated in cross-correlations above the threshold are collected and labelled as "multiplets".

Application to a fatigue crack growth test

This method has been applied to several acoustic emission dataset recorded during fatigue crack growth tests. Here, results from a load-imposed fatigue test ( R= σ min σ max =0.1 ) performed on compact tension (CT) specimen of a 5083 aluminium alloy are presented.

Two nano30 sensors, from PAC, were coupled on the specimen surface as depicted on fig. 4. The PCI2 system triggered recording waveforms sampled at 5 MHz when the signal exceeds the threshold of 43 dB. After the test, no filtering process was applied. Cross-correlation matrices and distributions were directly computed on the waveforms, independently for the both channels.

We selected a part of 1000 consecutive waveforms containing 2 multiplets for this study.

fig. 4 Geometry of CT specimen and sensors locations

Selection of the cross-correlation window

As explained in part 2.1, a proper window has to be set for the cross-correlation between two AE waveforms. In order to select the starting and ending point of time series, respectively n 0 and N, cross-correlation distributions are computed for different window lengths.

The pre-trigger of 20 µs is not taken into account in the calculation : n 0 is fixed, while N takes values to select [20, 40, 60, 80, 100, 120, 160] µs of the waveforms. Fig. 5 show an example of a waveform from the dataset and the different windows used for cross-correlations. The overall shape of the expected distribution presented in section 2.3 for all 7 windows can be found in Figure 6. By extracting signals involved in cross-correlations above 0.8, we clustered, as predicted, multiplets signals into one cluster.

But if the window is to short or too long, some multiplets signals are miss-classified. Extending the length of the cross-correlation window shifts the low cross-correlation peak towards 0 and modifies the large cross-correlation peak level. A cross-correlation performed on too short windows raises the risk to miss the right alignement resulting in an imprecise similarity measure, e.g. distributions computed on windows of 20, 40 and 60 µs whose maximums do not exceed 25 000 cross-correlations.. On the opposite, a too long window contains a large part of the waveforms codas which are less correlated and are characterized by a lower SNR. Hence, the distribution computed on windows of 160 µs has a lower maximum than distributions computed on windows of 80 and 120 µs. That is why the choice of the window is essentially ruled by the level of the large cross-correlation peak : the best clustering performance is obtain by maximizing the number of cross-correlation in the highest cross-correlation peak. Here, the window has to be set between [80,120] µs.

In addition, one observes that all distributions exhibit three peaks. One in low cross-correlation range and two in the large cross-correlation range. The second peak in the large cross-correlation range between 0.92 and 0.97 in fig. 6, not presented in part 2.3, contains cross-correlations performed between waveforms belonging to different multiplets emitted by the same type of sources. Selecting a threshold between these two close peaks allows to separate different multiplets thanks to an algorithm (not detailed here) by adding information of the signals emission time.

Noise effect

The noise robustness of our method is evaluated by adding white Gaussian noise to all waveforms. Distributions of cross-correlations are computed for different levels of SNR and are represented on fig. 7.

Fig. 7 Cross-correlation distributions for different SNR

From SNR = 100 to SNR = 0.1, the global shape of the distribution including high and low cross-correlations is preserved and shifted towards lower cross-correlation values, as expected. While the shape of the peak of low cross-correlations remains essentially unchanged, the two peaks of high cross-correlations are widening and merge. The method then is still applicable for low SNR, down to 0.1, by moving the cross-correlation threshold to separate highly similar waveforms, belonging to multiplets, from uncorrelated ones . Nevertheless, below SNR = 1 many errors are reported in the classification and below SNR = 5, it is not possible to distinguish the two peaks of high cross-correlations.

6.

Conclusion

This paper presents a method to automatically cluster groups of highly similar signals by computing the distribution of the cross-correlation matrix. This matrix is obtained by extracting the maximum of the cross-correlation function of each pair of waveforms from a fatigue test. The selection of the cross-correlation window is discussed and optimized by applying this method on several window lengths. And we show that this method is efficient to extract mulitplets, acoustic signature of fatigue crack growth, even if waveforms are polluted by noise.

Fig. 1 :

 1 Fig. 1 : Acoustic multiplets as a signature of fatigue crack growth (from [4]) -a) Waveforms of acoustic signal belonging to a multiplet during fatigue of Aluminum at 0.1 Hz, and location of the corresponding AE signals on the specimen ; b) Striations on a post-mortem SEM image of the crack face of the sample ; c) Acoustic activity during a fatigue test on aluminum : stress vs number of cycles. Magenta, cyan and orange clusters correspond to typical examples of multiplets ; d) Enlargement on a specific multiplet and e) on some loading cycles.

Fig. 5

 5 Fig. 5 Example of a waveform and the different cross-correlation windows. Each window starts at t = 0 µs


Bastian Bergfeld
email: bastian.bergfeld@slf.ch

Alec Van Herwijnen

Grégoire Bobillier

Eric Larose

Ludovic Moreau

Bertil Trottet

Johan Gaume

Janic Cathomen

Jürg Dual

Jürg Schweizer

Crack propagation speeds in weak snowpack layers

For the release of a slab avalanche, crack propagation within a weak snowpack layer below a cohesive snow slab is required.

As crack speed measurements can give insight into the underlying processes, we analysed three crack propagation events that occurred in similar snowpacks and covered all scales relevant for avalanche release. For the largest scale, up to 400 m, we estimated crack speed from an avalanche movie, for scales between 5 and 25 meters, we used accelerometers placed on the snow surface, and for scales below 5 meters, we performed a Propagation Saw Test. The mean crack speeds ranged from 36 ± 6 to 49 ± 5 m s -1 , and did not exhibit scale dependence. Using the Discrete Element Method and the Material Point Method, we reproduced the measured crack speeds reasonably well, in particular the terminal crack speed observed at smaller scales. Finally, we used a Finite Element model to assess the speed of different elastic waves in a layered snowpack. Results suggest that the observed cracks propagated as mixed mode closing cracks and that the flexural wave of the slab is responsible for the energy transfer to the crack tip.

INTRODUCTION

Dry-snow slab avalanche release is a complex critical phenomenon covering a wide range of scales (Schweizer and others, 2003). Across these scales, a sequence of fracture processes results in avalanche release. Dynamic crack propagation through a weak snowpack layer at the slope scale is one of these fundamental processes. Dynamic crack propagation has to be selfsustaining, so that the crack propagates below the slab across the slope; the slab will detach and an avalanche is released, if the slope angle is steep enough that gravitational pull overcomes frictional resistance. In other words, once dynamic crack propagation has started, slope topography and/or crack arrest will determine avalanche release size. Avalanche size is key in estimating the avalanche danger level [START_REF] Meister | Country-wide avalanche warning in Switzerland[END_REF], Statham and others, 2018, Techel and others, 2020). Reasons for crack arrest can be manifold. From a theoretical point of view, crack arrest occurs when the instantaneous dynamic energy release rate G falls below the dynamic fracture energy 𝛤 (e.g., Freund, 1990, p. 394): 𝐺(𝑟, 𝑟, ̇𝑡, loading, configuration, moduli, … ) < 𝛤,

1 2
Here G accounts for effects of loading, geometry, bulk material properties, instantaneous crack length r, time t and depends explicitly on the crack speed c = 𝑟̇ [START_REF] Freund | Dynamic Fracture Mechanics[END_REF]. The dynamic fracture energy 𝛤, a property characterizing the resistance to crack growth, may also depend on crack speed. Analytical solutions providing a dynamic energy release rate for crack propagation in weak snowpack layers are not available, and many of the required material properties remain unknown. What remains, from a theoretical point of view, is the suggestion that crack speed is an important factor which might be related to crack propagation distance.

Deriving an analytical solution for dynamic crack propagation in snow is a far cry, especially since the crack propagation mode is still debated. The crack propagation phase can be seen as a pure shear fracture process [START_REF] Mcclung | Shear fracture precipitated by strain softening as a mechanism of dry slab avalanche release[END_REF][START_REF] Mcclung | Fracture mechanical models of dry slab avalanche release[END_REF][START_REF] Mcclung | Application of dynamic fracture mechanics to dry snow slab avalanche release[END_REF], Bazant and others, 2003), while other studies suggested a mixed-mode anti-crack propagation process (Heierli and others, 2008, Gaume and others, 2018[START_REF] Mulak | Numerical investigation of the mixed-mode failure of snow[END_REF], Rosendahl & Weissgraeber, 2020b). So far, the debate is still open since there is no conclusive experimental evidence for either approach. Moreover, the fracture mode may depend on slope angle (Gaume andothers, 2019, Trottet andothers, 2021), the orientation of the crack relative to the slope, or possibly snowpack properties. The common fracture modes I -III are well defined for cracks in homogeneous materials, yet snow is a layered material. In snow, cracks generally propagate in thin weak layers, yet the energy required to create new crack surfaces (in the weak layer) stems from the slab layers above the weak layer. While equation 1 defines the balance between energy release (mainly in the slab, source of energy) and energy consumption (mainly in the weak layer, sink of energy), the energy has to be transferred from the source to the sink. Hence, there is an energy flux from the stress-strain field in the slab to the moving crack tip in the weak layer [START_REF] Broberg | Cracks and Fracture[END_REF].

The relative magnitude of the strain components in the slab defines the fracture mode and gives rise to the associated wave which transfers the energy through the slab to the crack tip. A crack propagation mode is therefore linked to an elastic wave, and the speed of the crack tip is therefore bounded by the speed of the associated type of wave [START_REF] Broberg | How fast can a crack go?[END_REF]. Crack speed measurements are therefore crucial, as these may provide an indication of the fracture mode.

Prior to 2000, no field measurements on crack propagation in snow were available. A few observations on whumpfs and firn quakes were published, and these phenomena were described as collapsing waves which can travel over large distances with speeds ranging from 6 m s -1 to slightly slower than the speed of sound in air [START_REF] Truman | Wave propagation in snow[END_REF][START_REF] Denhartog | Firn quake: A rare and poorly explained phenomenon[END_REF][START_REF] Benson | Stratigraphic studies in the snow and firn of the Greenland ice sheet[END_REF].

The first crack speed measurement in snow was reported by Johnson and others (2004), who used seismic sensors to measure the snow surface displacement during an artificially triggered, propagating crack in flat terrain, a so-called whumpf. They needed three days and numerous field experiments to obtain a single speed estimate, showing how challenging such in-situ measurements are. The mean propagation speed was 20  2 m s -1 over a propagation distance of about 8 m. Since then, numerous studies estimated crack speeds using high-speed photography of Propagation Saw Tests (PST, Fig. 1 upper left), a fracture mechanical field experiment for snow [START_REF] Sigrist | Critical energy release rates of weak snowpack layers determined in field experiments[END_REF][START_REF] Gauthier | Towards a field test for fracture propagation propensity in weak snowpack layers[END_REF]. Analysing displacements with Particle Tracking Velocimetry (PTV) provided new insight into weak layer fracture and crack propagation [START_REF] Van Herwijnen | High-speed photography of fractures in weak snowpack layers[END_REF], van Herwijnen and others, 2010, Schweizer and others, 2011, van Herwijnen and others, 2016b). Reported crack speeds range from 10 to 50 m s -1 (van Herwijnen & Birkeland, 2014[START_REF] Van Herwijnen | Fractures in weak snowpack layers in relation to slab avalanche release[END_REF], van Herwijnen and others, 2016b). Recently, high-resolution crack speed estimates were obtained with a digital image correlation (DIC) method, highlighting variations in crack speed along a PST and pronounced edge effects (Bergfeld andothers, 2018, Bergfeld andothers, 2021).

At scales larger than the typically 2 to 3 m long PST experiments, van Herwijnen and Schweizer (2011) derived a speed of 42  4 m s -1 over a distance of 60 m using seismic sensors deployed in an avalanche starting zone. To estimate crack speed, they determined the time between the first arrival of the P-wave and the avalanche release signal and assumed the propagation distance to be the width of the avalanche. Based on videos of avalanches, Hamre and others (2014) reported several widely varying estimates ranging from 18 m s -1 to 428 m s -1 . For the uncertainty of their speed estimates, they assumed an overall uncertainty on travel distance and travel time but distinct parameters (e.g. image resolution, propagation distance projected onto the terrain) of the individual experiments are not considered. Moreover, uncertainty is reported for a few speed estimates only. Hence, both studies suffer from methodological shortcomings, affecting the reliability of these crack speed estimates.

Overall, the most abundant and reliable data come primarily from PST measurements, raising the question whether the smallscale, rather one-dimensional PST experiment is representative of the two-dimensional, large-scale dynamic crack propagation that precedes slab avalanche release.

Performing experiments in avalanche terrain is challenging and often not possible due to safety concerns, so it is a good idea to apply numerical models to simulate this process. Recently, models based on the Discrete Element Method (Gaume and others, 2015, DEM;Gaume and others, 2017, Bobillier and others, 2020, Bobillier and others, 2021) and the Material Point Method (MPM; Gaume andothers, 2018, Gaume andothers, 2019) were developed to investigate crack propagation in snow.

Both methods were validated against PST experiments, showing their ability to reproduce crack propagation processes at the scale of a PST. While these numerical studies also provide new insight into crack propagation at scales beyond the PST scale, validation against slope scale experiments as well as an inter-comparison of both models using the same field experiment is missing.

As outlined before, the speed of cracks propagating in different fracture modes is limited by the speed of different elastic waves [START_REF] Broberg | How fast can a crack go?[END_REF] which can be computed from elastic properties and density of the snow. However, the layered snowpack, necessary for crack propagation acts as a waveguide that alters the speed and type of elastic waves travelling within. The Finite Element Method (FEM) can be used to study characteristics of guided wave propagation (Moreau and others, 2006[START_REF] Moreau | The use of an orthogonality relation for reducing the size of finite element models for 3D guided waves scattering problems[END_REF][START_REF] Castaings | Finite element model for waves guided along solid systems of arbitrary section coupled to infinite solid media[END_REF]. Such analysis is of great importance for non-destructive testing or health monitoring of engineering structures or the interpretation of seismic signals (e.g. [START_REF] Zhu | A unified approach for the structural health monitoring of waveguides[END_REF]Rizzo, 2012, Moreau andothers, 2020).

Hence, it seems also suitable for investigating wave modes in a layered snowpack.

Given the relevance of crack speed, our aim is to measure crack speed at a range of scales. To this end, we analyse three crack propagation events covering propagation distances from less than a meter up to more than 400 m. On the snowpack scale, we performed a PST experiment, where a 30 cm wide column is isolated and an artificial cut is introduced within a weak snow layer until crack propagation starts when the critical cut length rc is reached. Close to the PST, we triggered a whumpf, which is essentially a slab avalanche on flat terrain. At the whumpf, we measured crack speed with accelerometers placed on the snow surface. At the largest scale, we analysed a movie of an explosive-triggered avalanche that released at approximately the same time we performed the whumpf and the PST experiments but around 190 km away. Where possible, we reproduced these field experiments with numerical models based on DEM and MPM. To assess the upper limits of purely elastic guided waves for different propagation modes and to investigate the influence of slab layering on wave speeds, we used a FE model.

METHODS

The propagation events we analysed to derive crack speed consisted of a PST, a whumpf and an artificially triggered avalanche, covering distances from less than one meter to more than 400 meters and were recorded on 15 January 2019. For each of the events, we subsequently describe the analysis methods and estimate the snowpack properties since these serve as input for the numerical models (Fig. 1).

Field measurements

At each experimental site the snowpack was characterized with a manual snow profile following Fierz and others (2009) (Fig. 2). At the avalanche site, the manually observed profile was taken two days after the event and 200 m north of the avalanche in a slope of similar elevation, aspect and slope angle as the avalanche slope. For the PST, the crack propagated in a weak layer consisting of surface hoar (g, 10-15 mm) buried at a depth of 1.09 m (Fig. 2a, red dashed line) below a slab mainly consisting of decomposing and fragmented precipitation particles (c) and rounded grains (x). The crack at the whumpf site propagated at the top and near the avalanche site at the bottom of the weak layer that consisted in both cases of faceted crystals (e, 1-2 mm, Fig. 2b,c; red dashed lines). Slab thickness was 0.91 m and 0.88 m for the whumpf and avalanche, respectively.

To investigate variations in snow properties, we also performed snow micro-penetrometer (SMP) measurements at the PST and whumpf site. Generally, the heterogeneity within the PST was small (Fig. 3a) and the SMP measurements were in good agreement with the manual profile. A rather similar snowpack structure was observed at the whumpf site, where we took a SMP measurement at each sensor location; the main differences were a crust just above the weak layer at the whumpf site and the type of weak layer (surface hoar crystals vs. faceted crystals). The slab layering, identified in the manual profile (Fig. 2b), was manually transferred to the penetration resistance profiles to track snowpack changes along the crack propagation path (Fig. 3b, background colours).

Propagation Saw Test

We performed a 5.35 m long PST experiment on a flat and uniform site close to Davos, Switzerland. The experiment resulted in END and had a critical cut length rc = 38 cm. The exposed side wall of the PST was speckled with black ink (Indian Ink, Lefranc & Bourgeois) before the weak layer was cut with a 2 mm thick snow saw. A high-speed camera (Phantom, VEO710) recorded the speckled wall with an image resolution of 1280 pixel by 352 pixel and a rate of 7000 frames per second (lens aperture: f 2.8, focal length: 24 mm). Camera distortion correction and DIC analysis were done as described in Bergfeld and others (2021).

The DIC analysis (Subset size: 12 pixel, step size: 3 pixel, smoothing window: 181 frames) of the PST experiment provided us with displacement fields with time. In a post-processing step, we applied a threshold value of 0.05 mm (≈ five times the standard deviation of the noise in z-displacement before crack propagation) to the z-displacement to locate the crack tip at each time step (video frame) to derive crack speed, as the slope of linear fits in beam sections, along the PST (beam section width: 50 cm, step size: 5.3 cm). A detailed description of the applied methodology, including uncertainty estimates, is given in

Snowpack properties

For the PST experiment, slab, weak layer, and substratum thickness were measured in the manual profile (Fig. 2a). Density of the slab was measured in the field using a 100 cm 3 cylindrical density cutter (38 mm diameter) and a vertical resolution of 4 cm (Fig. 2a,b; black solid lines). Slab density is given in Table 1 as the mean of all layers attributed to the slab. Uncertainty is assumed to be 5 % (Proksch and others, 2016). Since density was not measured down to the ground, mean density of the substratum was estimated using a parametrization on hand hardness index and grain type suggested by [START_REF] Geldsetzer | Estimating dry snow density from grain form and hand hardness[END_REF].

The collapse height Δh and the effective elastic moduli of the slab Esl and weak layer Ewl were estimated from digital image correlation analysis. Collapse height was taken as the mean of the settlement of the slab measured after crack propagation. For the elastic moduli, 220 displacement fields with increasing crack length prior to crack propagation were used to compare to displacement fields predicted by the mechanical model established by Rosendahl and Weissgraeber (2020a). The optimal set of Esl and Ewl was estimated by minimizing the residual, with Esl and Ewl as free fitting parameters. Uncertainty is estimated as the standard deviation of the scatter over the different cut lengths. For a more detailed description see Bergfeld and others (2021). Since we did not have direct elastic measures of the substratum, its mean density was used to estimate the elastic modulus. To do so, the density parametrization of Gerling and others (2017) was scaled to fit with our PST measurement:

𝑬 = 𝟐. 𝟓𝟒 × 𝟏𝟎 𝟒 𝐤𝐠 𝐦 -𝟑 𝑨 ( 𝝆 𝝆 𝒊𝒄𝒆 ) 𝟒.𝟔 2
with the scaling parameter 𝐴 = 𝐸 sl PST 𝜌 sl PST ⁄ . In other words, we forced the density parametrization to fit our elastic modulus (𝐸 sl PST )density (𝜌 sl PST ) datapoint, determined for the slab in the PST experiment. Finally, we used a fixed Poisson's ratio of 0.25 (±0.05) for the slab, weak layer and substratum (Reuter and others, 2015).

Whumpf

Crack propagation also occurs on slopes not sufficiently steep to release an avalanche, a phenomenon called a whumpf (Johnson and others, 2004). To measure crack speeds in whumpfs, we used custom designed wireless time-synchronized accelerometers (ADXL362, cut-off frequency 100 Hz), with a resolution of 1 mG and a sampling rate of 400 Hz, placed on the snow surface. These sensors measure the acceleration of the slab when a crack in the weak layer passes by (Fig. 4). Though triggering whumpfs is difficult and unpredictable, we performed a successful experiment with seven sensors placed over a distance of 25 m on the same day we performed the PST experiment, and at a distance of 400 meters from the PST site. We aligned the sensors in a straight line and triggered a whumpf by jumping on the snow surface with skis. After the whumpf,

we measured the distances between the trigger point and each sensor using a tape measure. With the assumption of circular crack propagation, the crack propagation distance between two sensors Δx was taken as the difference in distance to the trigger point (Fig. 4b). To determine the onset of the acceleration in the i th sensor 𝑡 𝑖 𝑜𝑛𝑠𝑒𝑡 (dots in Fig. 4b) we used the Akaike Information criterion AIC (Kurz and others, 2005). The time difference Δt = 𝑡 𝑖+1 𝑜𝑛𝑠𝑒𝑡 -𝑡 𝑖 𝑜𝑛𝑠𝑒𝑡 between two sensors was then used to compute crack propagation speed as ci,i+1 = Δ𝑥 Δ𝑡 . Uncertainty was assessed by Gaussian error propagation and the contributions from the distance measurement (ux = 0.1 m) and the onset determination. For the latter, a systematic uncertainty was assigned by visually inspecting each AIC automatic onset time on the waveform and comparing it to a manual onset estimate.

The frequency content of the acceleration signals was computed with consecutive Fourier transforms (100 samples each, overlap 90%) using SciPy 1.5.0 (Virtanen and others, 2020) (Fig. 5). Most of the energy was below 50 Hz, and the median central frequency over all seven recorded whumpf signals was 7.6 Hz.

Snowpack properties

Slab, weak layer, and substratum thickness were measured in the manual profile (Fig. 2b) and shown in Table 1. For density we used a cylindrical density cutter and sampled the slab layers every 4 cm. There were no density measurements for the weak layer and substratum. We thus again used the approach of [START_REF] Geldsetzer | Estimating dry snow density from grain form and hand hardness[END_REF] to estimate its mean density. Since we cannot measure elasticity in the field, mean density of the substratum and slab was used to estimate the elastic modulus using Equation 2. Values of elastic moduli of weak layers are very scarce in literature, and we had no possibility to measure it directly. Instead, we assumed the elastic properties to be the same as estimated in the PST experiment. The displacement of the slab was assessed by integrating the acceleration signal twice and forcing the velocity signal (after first integration) to be zero after collapse (Fig. 5a, orange line). The collapse amplitude of the whumpf Δh = 2.2 ± 0.6 mm was taken as the mean of the end displacements of sensors three to seven, since those signals were not influenced by the triggering.

Avalanche

We analysed a video of an explosive-triggered avalanche that released the day we performed the other two experiments. The avalanche released on a steep slope (~40°) above the ski resort of Grimentz, Switzerland (WGS84: 46.1752° N, 7.5200° E), located about 190 km from the field site in Davos. The avalanche was filmed at 30 frames per second and an image resolution of 1920  1080 pixel 2 . For georeferencing surface cracks and avalanche outlines, we analysed selected movie frames with the Monoplotting tool (Bozzini and others, 2012) to obtain georeferenced vector data by drawing directly on the image (Fig. 6 left and middle). Transferring the vector data to a GIS-system, we then estimated the crack propagation path by taking the line of sight between the trigger point (location of the bombing) and the location where the surface crack opened, whereby visible rocks were bypassed. Crack propagation distance was then computed by projecting the path onto the terrain and computing the distance from the trigger point to the surface crack. The time the crack needed to propagate this distance was estimated from the movie as the lag time between the detonation frame and the frame in which the associated surface crack became visible. Uncertainty of the crack propagation speed uc was computed with Gaussian error propagation and originated from an uncertainty in distance ux = 0.05x (5% of the distance, and at least 5 m) and time ut = 5 frames = 0.17 s. For each crack propagation path, we further derived slope angle and crack orientation relative to the slope aspect from a digital elevation model (5-meter resolution, source: Swiss Federal Office of Topography swisstopo, Berne, Switzerland). A mean slope angle of every crack path was estimated by extracting the slope angle every three meters along the path. Every 3 meters along the path, we also took the azimuths of the crack orientation and the slope aspect. The difference between these two angles was normalized to show if the crack propagated down-, cross-or up-slope.

Snowpack properties

Slab, weak layer, and substratum thickness are shown in Table 1 and were taken from the manual profile (Fig. 2c). Layer densities were estimated using the approach presented by [START_REF] Geldsetzer | Estimating dry snow density from grain form and hand hardness[END_REF]. Applying equation 2, we estimated the elastic modulus using mean density of the substratum and slab. Concerning the weak layer, the elastic modulus was again assumed the same as in the PST experiment, and the collapse amplitude Δh = 4 mm was taken as the median collapse amplitude from 192 experiments reported by van Herwijnen and others (2016a). (1979), is composed of a large number of discrete connected and interacting particles. The PST system was generated using commercial DEM software; PFC3D (v5)

(http://www.itascacg.com). For a detailed description of the DEM model for snow, and in particular how particle (and particle contact) parameters relate to macroscopic snow parameters, refer to Bobillier and others (2020). Recently, Bobillier and others (2021) modelled the PST design and showed that crack propagation in weak snow layers can be reproduced with their DEM model (Fig. 7a). They simulated a three-dimensional PST consisting of three layers: a basal layer, a weak layer and a uniform slab layer on top. Where available, we used the snow parameters measured in the field (Table 1), otherwise we used the values proposed by Bobillier and others (2021) (see Appendix A), obtained by matching simulation and experimental results. To derive crack propagation speed, Bobillier and others (2021) compared different methods, based on weak layer stress, weak layer bonding damage and slab displacements. They concluded that all methods provide comparable results. In accordance with our field measurement, we used their method based on slab displacements to estimate crack propagation speed in the DEM simulation of the PST. Due to computational limitations the PST is the only field experiment which could be reproduced with the DEM.

Material Point Method

The MPM method is a continuum and hybrid Lagrangian-Eulerian numerical technique (Sulsky and others, 1994), making it particularly well-suited to handle processes involving large deformations and fractures. Snow constitutive models based on critical state soil mechanics [START_REF] Roscoe | On the generalized stress-strain behaviour of 'wet' clay[END_REF] were developed to simulate the mechanical behaviour of snow slab and weak snow layers [START_REF] Stomakhin | A material point method for snow simulation[END_REF]others, 2013, Gaume andothers, 2018). Gaume and others (2018) reproduced results from different PST experiments, and then applied their model to simulate slab avalanche release on an ideal slope. Here, the numerical setup consisted of a PST with a uniform slab and a weak layer. Particles at the bottom of the weak layer were fixed in position (Dirichlet boundary condition). Details about the constitutive models can be found in Gaume and others (2018) and Trottet and others (2021). Where available, we used the snow parameters measured in the field (Table 1), otherwise we used the values which were proposed for snow by Gaume and others (2018) (see Appendix B), obtained by matching simulation and experimental results. To initiate crack propagation in the simulation, an initial crack of length r was generated in the weak layer by setting the elastic modulus of the associated particles to zero. Quantities of interest, such as the position x or the volumetric plastic strain, were stored for each particle. The crack length was then tracked by defining the crack tip as the location of the furthest plasticized particle, and crack speed c was directly obtained from the temporal evolution of the crack tip position.

Elastic limits and crack speed estimates

Crack speed estimates

The speed of propagating cracks in snow was estimated by [START_REF] Heierli | Solitary fracture waves in metastable snow stratifications[END_REF] , g the gravitational acceleration, h the collapse height, 𝜌 the mean slab density, H the slab thickness and ν the Poisson's ratio (0.25 ±0.05). Basically, Heierli (2005) followed Johnson and others ( 2004) who reported a measurement of crack propagation speed in flat terrain and suggested that their measurement should be explained as a flexural wave in the slab. [START_REF] Mcclung | Approximate estimates of fracture speeds for dry slab avalanches[END_REF] proposed an alternative explanation for the measurement of Johnson and others (2004) and estimated empirical terminal crack speeds to be in the range 𝑐 sc = (0.7 -0.9) 𝑐 s where 𝑐 s is the shear wave speed: 𝑐 s = √ 𝐸 𝑠𝑙 2(1+ν)𝜌

. Uncertainties for the crack speed estimates are consequently propagated from the uncertainties of the required snowpack parameters (Table 1).

Elastic wave speeds

Calculating the energy flux into the crack tip region leads to upper crack speed limits [START_REF] Broberg | How fast can a crack go?[END_REF]. For mode I and mode II cracks, the Rayleigh wave speed 𝑐 𝑟 = √ 0.87 • 1.12ν 1+ν 𝑐 𝑠 is the upper limit (Chi [START_REF] Vinh | An approach for obtaining approximate formulas for the Rayleigh wave velocity[END_REF]. If the forbidden subsonic super-Rayleigh region is bypassed in mode II cracking, then the upper limit is the longitudinal wave velocity [START_REF] Broberg | How fast can a crack go?[END_REF]. A mode III crack is limited by the shear wave speed 𝑐 𝑠 . Again, uncertainties are propagated from input parameters. However, the above expressions are valid for waves travelling in an elastic full or half space, whereas the natural snowpack is a layered medium and acts as waveguide.

𝑐 𝑙 = √ 𝐸 𝑠𝑙 𝜌

Finite Element Method

A snowpack associated with slab avalanche release consists of at least three layers: slab, weak layer, and substratum (Fig. 8d).

Due to the limited thickness of these layers, the entire system acts as a waveguide, affecting wave modes, propagation speed and damping. To investigate which wave modes travel at what speed along the length of the layered slab, FE simulations of wave propagation in a PST-like setup were made with the commercially available software COMSOL Multiphysics. The extent of the model in the x direction was set to 300 m, in order to avoid interferences between the direct wave and unwanted boundary reflections, which would cause resonances of the entire structure. The domain was meshed with squared elements of one meter in size. The source was defined as a vertical force located at the free surface of the slab, at x = 0 m (Fig. 8a). The source was modelled as 2-cycles gaussian tone burst with central frequency of 50 Hz. From the energy of the wavefield, different modes can be distinguished (Fig. 8c, e.g. red and blue boxes). In a homogeneous waveguide with only one layer, instead of three, the wavefield consists of a superposition of the so-called Lamb modes. When the frequency is higher than the cut-off frequency of the higher order modes, many modes can propagate at the same time, and the energy of the wavefield splits between the modes, with a distribution that depends on the frequency. When the frequency remains under the cut-off frequencies, only the two fundamental modes propagate in the waveguide. At very low frequencies, the asymptotic behaviour of these two modes can be approximated by a flexural and a longitudinal wave.

In the present case of a waveguide with three layers, however, each layer may act as a waveguide of its own that is coupled with the other two layers. As such, there is co-existence of more than two modes, even from an asymptotic point of view.

There is energy exchange between the modes, with a distribution that depends not only on the frequency, but also on the variations of acoustic impedance between the layers. This can be seen for example in Fig. 8c, in the region marked with a red box. The frequency-wavenumber branch in this box looks like a flexural mode with a discontinuity around 5 Hz, but the energy of the wavefield is actually split between two co-existing flexural-like modes. The same goes for co-existing longitudinal-like modes, which can be seen in the blue box. The figure also shows branches of higher order modes, with a cut-off frequency around 50 Hz.

For simplicity, we will associate in the following the branches in the red and blue boxes, in Fig. 8c, as one flexural and one longitudinal mode, respectively. This approximation can be justified by the motion of these branches, which is flexural and longitudinal. The nonlinearity between the wavenumber and the frequency (Fig. 8c) indicates that the modes are dispersive.

Hence different frequencies travel at different velocities. A consequence of this dispersion is that a wave packet centered around a given frequency travels with a group velocity defined by 𝐶 𝑔 = 𝜕𝜔 𝜕𝑘 , where  is the angular frequency and k is the wavenumber. To calculate the group velocity from the frequency-wavenumber spectrum, we proceed in two steps: For each frequency, 𝑓 𝑛 , we find the corresponding wavenumber, 𝑘 𝑛 , by extracting the value of k where the spectrum amplitude is maximum. The group velocity, 𝐶 𝑔 (𝑓 𝑛 ) around the n th point (𝑓 𝑛 , 𝑘 𝑛 ) in the spectrum, is finally approximated numerically such that 𝐶 𝑔 (𝑓 𝑛 ) = 2𝜋 In most models, also in the DEM and MPM models, the slab is treated as a uniform layer with an isotropic effective elastic modulus and a mean density. In nature, the slab generally has a pronounced layering. In the density profile (Fig. 8b) we visually identified 6 sections where density approximately changed linearly with depth (black line in Fig. 8b), and converted those into a section-wise elastic modulus profile using Equation 2(red line in Fig. 8b).

To highlight the influence of slab layering on wave speed, we performed two simulations, one with a uniform slab and one with a layered slab. In addition, we investigated the influence of weak layer properties on crack speed by varying the elastic modulus of the weak layer (0.4 MPa, 0.6 MPa, 1 MPa). In total, we performed six FEM simulations.

RESULTS

We derived crack propagation speeds from three crack propagation events recorded on the same day. These experimental values were compared to results of numerical models established to reproduce crack propagation in weak snow layers (DEM, MPM) or to assess elastic wave speeds in the slab (FEM).

Field measurements

In the PST, crack propagation started with an initial speed of 17 m s -1 in the first meter and increased thereafter (Fig. 9a, blue line). Crack speed peaked at x = 4.5 m with 58 m s -1 , and overall the mean speed was 39±13 m s -1 . In the whumpf experiment, crack speed was computed for different sensor combinations. The two sensors closest to the trigger point could not be used, as the signal also contained noise from jumping on the snow (Fig. 4b, blue and orange line). Also, we did not compute the speed between sensors 6 and 7, as these were at a very similar distance to the trigger point (Fig. 4a, brown and violet sensors). In total, we thus had 9 pairs to calculate crack speed (Fig. 9a, orange dots), with a mean speed of 49±5 m s -1 . In the avalanche movie, we analysed 14 crack paths, resulting in crack propagation speeds between 23 and 44 m s -1 (mean 36±6 m s -1 ) covering distances from 26 to 440 meters (Fig. 9a, green squares). Crack propagation paths were located on slopes with mean slope angles between 32° and 42° and there was no link between slope angle and crack speed (Pearson r = -0.51, p = 0.05; Fig. 9b).

Since the explosive was placed near the ridge, there was only one upslope crack path, while all others were cross-slope or downslope. Again, there was no clear link between propagation path and crack speed (Pearson r = -0.13, p = 0.63; Fig. 9c).

Numerical modelling

Both the DEM and the MPM simulations reproduced the right order of magnitude as well as the overall trend in crack speed along the PST experiment. While values obtained with the DEM model underestimated crack speed compared to the experiment, the MPM model overestimated crack speed. Away from the PST beam ends (2 m < x < 4 m), where edge effects are small (Bergfeld and others, 2021), mean speeds were 44±6 m s -1 , 31.6±2.1 ms -1 and 63±3 m s -1 for the experiment, DEM and MPM, respectively.

At the whumpf scale, we used two different MPM geometries. First, we modelled a 3-dimensional PST-like beam configuration (width: 0.3 m, length: 25 m, labelled as "MPMbeam" in Fig. 11a). Second, we modelled a 3-dimensional flat area of 25 m by 25 m (MPM -areal in Fig. 11a). Model geometry had no influence on crack propagation speeds, as mean speeds were 54±2.2 m s -1 and 56±1.6 m s -1 for the areal and beam configuration, respectively, in good accordance with the field measurements (49±5 m s -1 ).

The avalanche was also modelled using two different MPM configurations. First, we simulated a cross-slope PST-like beam configuration (width: 0.3 m, length: 25 m, slope angle: 40°, labelled as "MPMbeam" in Fig. 11b). Second, we modelled a tilted area (40°) of 25 m by 25 m and computed the cross-slope crack propagation speed (MPM -areal in Fig. 11b). The beam configuration resulted in a fairly constant speed of 𝑐̅ = 48±4 m s -1 , in the same range as the speeds derived from the movie.

The slope simulation showed a steady increase in speed (Fig. 11b, purple solid line), without reaching a steady state within the 25 m of the model domain.

Elastic limits and crack speed estimates

Elastic wave speeds were in the range of 97±6 to 224±35 m s -1 (Table 2). As expected, they were generally much higher than the experimental crack propagation speeds. McClung (2005) estimated a realistic range for weak layer crack speeds to be csc = (0.7 -0.9) cs which would lead to csc ∈ [74, 137] m s -1 and therefore around double the experimentally estimated crack propagation speeds. The estimate for crack speed 𝑐 𝑓𝑤 as suggested by [START_REF] Heierli | Solitary fracture waves in metastable snow stratifications[END_REF], agreed well with the experimental crack propagation speeds 𝑐 𝑒𝑥𝑝 (Table 2). In our FEM analysis the PST-like beam configuration, consisting of slab, weak layer and substratum, acts as a waveguide with multiple layers. We determined the group velocity of the flexural and the longitudinal wave mode for different weak layer elastic moduli, as well as with and without slab layering. The flexural wave speed cf increased up to around 30 Hz (120 m s -1 ), and decreased thereafter (Fig. 12a). At higher frequencies (f > 35 Hz) the decrease was more pronounced for a layered slab (compare dashed and dotted lines in Fig. 12a). In the region of interest, at a dominant frequency of 7.5 Hz, group velocities for the flexural mode ranged from cf = 95 m s -1 to 105 m s -1 , indicating that the elastic modulus of the weak layer as well as slab layering had little influence on flexural wave speeds.

Overall, longitudinal wave speeds were much higher at low frequencies (f < 35 Hz) than flexural wave speeds. Weak layer elastic modulus did not influence the longitudinal wave mode. However, longitudinal wave speeds were somewhat sensitive to slab layering in the lower frequency range. Longitudinal wave speeds at 7.5 Hz were in the range of 370 m s -1 to 380 m s -1

for the layered slab and 440 m s -1 to 470 m s -1 for the uniform slab.

DISCUSSION

We presented crack speed measurements from three crack propagation events covering a wide range of scales and compared these to theoretical wave speed limits and results of numerical simulations. Overall, our measurements of crack speed from the PST, whumpf and avalanche slope provided similar values, although there were some differences between the different propagation events (Fig. 9). We did not observe a dependence of crack speed on propagation distance, suggesting that crack propagation in the whumpf experiment and the avalanche had a fairly constant speed around, 49 m s -1 and 36 m s -1 , respectively.

Determining crack speed always requires tracking the crack tip. Hence, any method to estimate the crack tip potentially affects crack speed estimates. In our flat field experiments, PST and whumpf, we estimated the location of the crack tip using vertical displacement and acceleration, respectively. The main source of energy for crack propagation stems from the vertical displacement of the slab, since horizontal displacement does not free any energy (e.g. gravitational potential energy or mechanical energy stored in the slab). Hence, the crack tip speed has to be in line/synchronous with the speed the vertical displacement travels along the slab. Therefore, crack speed estimates derived from vertical displacement and acceleration in flat field experiments provide accurate estimates, even if the tracked location does not necessarily coincide with the exact location of the crack tip. Hence, the methods to determine crack speed are robust in case of the PST and whumpf, the crack speed in the avalanche, however, could be greater than our methodology revealed. It remains an open question if slab fractures, visible at the snow surface, are good proxies for the location of the crack tip in the weak layer. Indeed, it is entirely possible that the crack in the weak layer has already propagated further before a crack in the slab becomes visible on the images. We considered three possible processes (Appendix C) that could result in a delay between the appearance of slab fractures and the passing of the crack tip in the weak layer, resulting in lower crack speed estimates. However, as the processes are presumably constant, the influence on crack speed decreases with increasing propagation distance. By accounting for these delays, we more accurately estimated the terminal crack speed (45 ± 5 m s -1 ), which was still close to the uncorrected estimate (36 ± 6 m s -1 ).

Since similar experimental crack speeds were observed throughout the three crack propagation events, it can be assumed that crack propagation processes observed in the PST were similar to those in the whumpf and the cross-slope propagation in the avalanche. The increase of crack speed in the first two meters of the PST beam highlights the need for long PST experiments, at least longer than the so-called touchdown distance, which is typically 2 to 3 m (Bair andothers, 2014, Bergfeld andothers, 2021). If this is fulfilled, our results suggest that crack propagation in a PST is representative for larger 2-dimensional crack propagation occurring in whumpfs and for the cross-slope direction in our avalanche movie. Long PST experiments may thus be well suited to investigate self-sustained crack propagation.

On the PST scale, the numerical simulations (MPM, DEM) reproduced the initial increase in crack speed well. In general, the DEM method underestimated crack speed compared to the experiment, while the MPM method resulted in higher values (Fig. 10). The discrepancy of the crack speeds between the MPM and DEM methods could stem from the different implementations how damping is realized in the models. In MPM, no explicit numerical damping is used at the level of the balance equations and the constitutive model is rate-independent. A small amount of damping is introduced during the interpolation between the grid and the particles (see Gaume and others (2018) for more details) but does not affect the presented results. In contrast, in DEM, a global damping coefficient of 0.6 [START_REF] Cundall | A discrete numerical model for granular assemblies[END_REF] was used in the Newton's equations.

The DEM method is computationally very costly (the PST took around one day on a 14 cores Intel Xeon CPU with 256Gb RAM) and is therefore currently limited to the PST scale. The MPM method, on the other hand, can also be used to perform simulations for geometries larger than the PST, provided certain requirements are fulfilled. Currently, the mesh size has to be the same for the modelled system containing slab and weak layer. Resolving processes in the weak layer therefore requires a mesh size smaller than the weak layer thickness. Reducing the mesh size was computationally too costly. Instead we modelled the weak layer with a thickness of 45 cm rather than the 1.5 cm observed in the manual profile (see Appendix B). Since the thickening of the weak layer affects the stiffness of the weak layer, we scaled the elastic modulus of the weak layer accordingly to keep the stiffness constant. Doing so, the measured crack propagation speed in the whumpf was fairly well reproduced, independent of the MPM model geometry (Fig. 11a). This MPM result further suggests that the PST-like configuration, which is rather 1-dimensional, is representative of 2-dimensional crack propagation in flat terrain as observed in whumpfs.

For the avalanche slope, we again tested two MPM model configurations, with very different outcomes. Here, the crack speed of the areal configuration increased with propagation distance and did not reach a terminal speed within the 25 m (Fig. 11b, purple solid line). Hence, in the areal configuration the MPM model predicted a different crack propagation speed behaviour than derived from the avalanche video. Probably, because the crack in the areal MPM simulation propagated in Mode III. In contrast, an increasing crack speed with propagation distance was not observed for the cross-slope, beam like MPM configuration (Fig. 11b, purple dashed line), which resulted in similar crack speeds as derived from the avalanche video. Most likely this is due the downside support of the beam, which we added to force the cracking in a pure flexural mode similar as in the whumpf and PST experiment, where mode III does not occur.

Therefore, we assume that the cross-slope crack propagation observed in our avalanche movie is not influenced by the slope angle, which potentially could change crack propagation into pure mode III. Possibly, consecutive slab fractures behind the crack tip hinder this transition into mode III. In the areal MPM simulation the slab is not allowed to fracture and consequently showed this transition, indicated by the increasing crack speed. Another factor which may have prevented the transition into mode III in our avalanche movie could be the natural supporting effect from the foot of the slope. Gaume and others (2019) used MPM to simulate crack propagation in a full slope, which reached into flatter areas as well, and did not report an increasing crack speed in the cross-slope direction. In contrast, our areal MPM configuration had a uniform slope angle and was truncated at the model boundaries, hence, downside support was not possible for geometrical reasons. The mismatch between crack speeds from the avalanche movie and our MPM simulation suggests that more research is needed to understand the influence of slope angle in the MPM model, and likewise, more field data are needed to see whether our single measurement is representative.

With the exception of the latter mismatch, we can say that the MPM and DEM derived crack speeds are in reasonable agreement with the measurements. However, given the large number of input parameters in the models, and that some of these parameters cannot (or only imprecisely) be measured in the field, a sensitivity study would be required to add some uncertainty to the model results, but this is beyond the scope of the present work.

Further, we compared our experimental crack speeds to the speed of body waves in infinite media and to two models predicting crack speed in relation to slab avalanche release (Table 2). Only the formulation suggested by [START_REF] Heierli | Solitary fracture waves in metastable snow stratifications[END_REF] agreed well with experimental values, while the crack speed range suggested by [START_REF] Mcclung | Approximate estimates of fracture speeds for dry slab avalanches[END_REF] and the body wave speeds were substantially higher. Of course, the latter is expected since typical maximum crack speeds of other materials are typically less than 60 % of the corresponding wave speed [START_REF] Broberg | How fast can a crack go?[END_REF].

As crack propagation preceding avalanche release occurs in a layered medium, we also performed FEM simulations to identify speeds of different wave modes as elastic limits for crack propagation in snow. Such a system can be seen as a waveguide which introduces geometric dispersion. In that case, the group velocity carries most of the energy and is dispersive (i.e. frequency dependent). A spectral analysis of the observed whumpf signal showed that most of the energy of the propagating crack was below 50 Hz (Fig. 5), with a median peak in the spectrum around 7.5 Hz. This is in good agreement with the spectral content of seismic signals generated by crack propagation before slab avalanche release [START_REF] Van Herwijnen | Seismic sensor array for monitoring an avalanche start zone: design, deployment and preliminary results[END_REF].

Evaluating the FEM simulated group velocities at 7.5 Hz, the measured crack speed was around 0.4 times the flexural wave speed, while it was only 0.1 times the longitudinal wave speed.

It has been known for a long time that cracks tend to accelerate to a constant (terminal) speed [START_REF] Schardin | Velocity effects in fracture[END_REF], which is lower than the elastic wave speed that corresponds to the crack propagation mode [START_REF] Broberg | How fast can a crack go?[END_REF], i.e. the elastic wave speed is an upper limit. This restriction originates from the energy flux into the dissipative region, namely the crack tip. The energy flux decreases with increasing crack speed, and it goes down to zero at the limiting speed [START_REF] Broberg | How fast can a crack go?[END_REF][START_REF] Freund | Dynamic Fracture Mechanics[END_REF]). In the avalanche case, the available energy stems from the work of gravity released by the slab, then the energy is transferred by the slab to the moving crack tip within the weak layer. Therefore, the crack propagation speed is limited by an elastic wave mode of the slab. Since the PST and whumpf experiments were performed on flat terrain, we assume that the energy available for crack propagation stems from flexural slab bending due to weak layer collapse. Since for self-sustained crack propagation the crack tip cannot propagate faster than the energy source, we suggest that the flexural wave speed is an upper bound for the crack speed on flat terrain. This also explains why the formulation of [START_REF] Heierli | Solitary fracture waves in metastable snow stratifications[END_REF] agreed very well with our measurements, as it is essentially a formulation for a flexural wave. In other words, if the fracture mode is a mixed mode anticrack (Heierli and others, 2008), its limiting speed is given by the flexural wave speed in the slab cf. In our experiments, the speeds were c ≈ 0.4 cf.

The similarity of the crack speed measured in the avalanche experiment suggests that crack propagation on the slope did not change the mode of propagation, whereas tilting the model geometry changed the mode in the MPM simulation. That the fracture mode did not change is also supported by the fact that we found neither a dependence of crack speed on crack orientation nor on slope angle (Fig. 9). However, these observations are based solely on the results from one single avalanche.

We do not demonstrate that other fracture modes do not exist in nature. In particular, for distant triggering through steep slopes, gravitational tensile stress might favour the crack propagation in mode II and III, and crack propagation at greater speeds than observed in our single measurement might occur as suggested by recent modelling results (Gaume andothers, 2019, Trottet andothers, 2021). Clearly, future measurements are needed to confirm this interpretation.

Our FEM analysis further revealed that slab layering did not overly influence the flexural wave speed and therefore, probably had little effect on the crack propagation speed. This is of particular interest since the slab is most often treated as homogeneous.

Another parameter hard to determine in the field is the weak layer elastic modulus, and its influence on the flexural wave speed was also very limited.

CONCLUSIONS

We presented three examples of crack speeds in weak snowpack layers determined from field experiments and observations, and compared these to results from numerical models and elastic wave speeds. Crack propagation speed is a highly relevant driver of dry-snow slab avalanche release size.

The three propagation events (PST, whumpf and avalanche) occurred on the same day in similar snowpacks and cover a range of distances from less than 1 m to more than 400 m, hence the scales relevant for slab avalanche release. The speeds ranged from 39 to 49 m s -1 , which are remarkably close in view of the natural variations in the different parameters, e.g. the slab thickness. The similar speeds obtained for different distances and geometries allow the assumption that crack speed measurements obtained from PST experiments may be representative of slope scale crack propagation, provided the columns are long enough, ideally greater than 5 meters. Crack speeds modelled with DEM and MPM were in the same range as the experimental values. The MPM simulation for a tilted area of 25 m  25 m revealed a steadily increasing crack speed, which was not observed in our single avalanche event. The comparisons also showed some discrepancies, demonstrating the need to better understand the link between model results, model input parameters and how these relate to snowpack properties.

Using a FEM model, we determined the dispersive speeds of different guided wave modes. The crack speeds we obtained corresponded to about 0.4 times the flexural wave speed at 7.5 Hz, the dominant frequency in signals measured during crack propagation with wireless accelerometers. The FEM results further suggest that slab layering and weak layer elastic properties may have little influence on wave speeds at this frequency.

Our measurements of crack speed at various scales are limited, yet they are consistent with the assumption that crack propagation in weak snowpack layers, at least in flat terrain, is linked to the flexural wave mode in the slab in the same way that a mode II crack in other materials is generally linked to elastic shear or longitudinal waves. Measurements on steep slopes and in different snowpacks, as well as comprehensive modelling approaches, are still needed for solid conclusions on the dominant fracture modes in dry-snow slab avalanche release.

Author contributions. JS, AH and BB designed the research before BB and AH conducted the field experiments. GB performed the DEM simulations, BT and JG made the MPM simulations and LM together with EL did the FEM analyses. JC helped to geo-reference the avalanche movie. JD and JS contributed to all parts of the project. The manuscript was written by BB with input from all co-authors.

Fig. 1 .

 1 Fig. 1. Overview of methods and scales used to investigate crack speed. Top row: The experimental methods to measure crack speed consisted of a PST, a whumpf and an artificially triggered avalanche, covering distances from less than one meter to more than 400 meters. Bottom row: The numerical models used to reproduce these experiments were based on the Discrete Element Method (DEM), the Material Point Method (MPM) and the Finite Element Method (FEM).

Fig. 2 .

 2 Fig. 2. Manual snow profile showing hand hardness index (width of bars), grain type (colours, see legend) and layer density (black line) for (a) PST, (b) whumpf and (c) avalanche site. Layer density was either measured with a density cutter (black solid line) or was estimated using grain type and hand hardness index (black dashed line). The depth of the crack in the weak snowpack layer is indicated with the red dashed lines.

Fig. 3 .

 3 Fig. 3. (a) Penetration resistance measured with the SMP along the PST experiment (blue lines plotted on an image of the speckled side wall of the PST). (b) Penetration resistance measured at each acceleration sensor (indicated by colour as shown in Fig. 4) after triggering the whumpf. The background colours indicate the grain type (top right) of the corresponding layer in the manual profile. Measurements are restricted to the slab to the top of the melt-freeze crust (red layer in Fig. 2b).

Fig. 4 .

 4 Fig. 4. (a) Arrangement of the accelerometer sensors at the whumpf site. (b) Measured downward acceleration with time. The time difference Δt of the onsets of acceleration (black dots) between sensor pairs in combination with their spacing Δx were used to compute crack speed.

Fig. 5 .

 5 Fig. 5. (a) Downward acceleration (blue line) and displacement obtained by double integration (orange line) with time for sensor 4 at the whumpf site (see Fig. 4). (b) Corresponding spectrogram showing the energy (colours) per frequency band with time.

Fig. 6 .

 6 Fig. 6. Schematic showing the work flow for the crack speed estimation based on the movie. (left) Selected frames from the movie were exported to analyse the position of tensile cracks. (middle) the Monoplotting tool was used for geo-referencing cracks in the frame. (right) These geodata were transferred to a map for estimating crack paths and the corresponding projected crack propagation distances (map source: Federal Office of Topography swisstopo).

Fig. 7 .

 7 Fig. 7. PST setup of the (a) DEM and (b) MPM modelling while sawing the PST. The colouring in the slab indicates the vertical displacement. The insets (15 cm x 15 cm) show a close-up around the saw. The DEM inset (a) shows slab and weak layer particles in blue and green, respectively. Around the saw particles are clipped and bonding damage in the weak layer is highlighted in red. In the MPM inset (b), the crack (red) started propagating and is already ahead of the snow saw.

  [START_REF] Mcclung | Approximate estimates of fracture speeds for dry slab avalanches[END_REF].[START_REF] Heierli | Solitary fracture waves in metastable snow stratifications[END_REF] proposed an analytical model, suggesting that the crack in the weak layer occurs in the form of a localized disturbance zone (crack tip to the trailing point where the slab rests on the substratum again) propagating as a flexural wave with constant speed 𝑐 𝑓𝑤 = √-

Fig. 8 .

 8 Fig. 8. Setup of the FEM simulations. (a) The propagation of guided waves in the slab of a PST-like arrangement was simulated by inducing a pulse at x = 0 m and using (b) an elastic modulus profile for the slab layer. (c) The dispersion relations of the different wave modes were used to compute their propagation speeds. (d) Side view of the PST-like FEM configuration consisting of the substratum, weak layer and slab.

Fig. 9 .

 9 Fig. 9. (a) Crack speed with propagation distance for the PST (blue), the whumpf (orange) and the avalanche (green). For the PST, the red area behind the blue curve indicates the uncertainty. For whumpf and avalanche, the uncertainty is given with red error bars. The grey horizontal lines behind the whumpf speeds indicate the distance range used to compute the respective data point. (b) Crack propagation speed for the avalanche with mean slope angle and (c) with crack orientation relative to the slope aspect. Each dot represents a speed estimate from a cracking path while the bars show the uncertainty.

Fig. 10 .

 10 Fig. 10. Crack speed estimates for the PST. The Discrete Element Method (DEM, brown line), the Material Point Method (MPM, purple dashed line) and the field experiment (blue) showed a strong increasing crack speed at the beginning of crack propagation before crack speed increased less in the centre part of the PST beam

Fig. 12 .

 12 Fig. 12. Group velocities of (a) the flexural wave mode and (b) the longitudinal wave mode with frequency. Simulations with a uniform effective elastic modulus for the slab are shown as dashed lines, while the crosses indicate a layered slab with varying elastic modulus as shown in Fig. 8b. Different values for the elastic modulus of the weak layer are depicted in black, red and blue.

Table 1 .

 1 Measured or estimated snowpack properties for the PST, whumpf and avalanche. The columns D, M and F indicate if the property is used as input for the DEM, MPM or FEM simulations, respectively.

	Experiment	PST	D M F	Whumpf	M	Avalanche	M
	Slab thickness (m)	1.09 (±0.03)	x	x	x	0.91 (±0.03)	x	0.88 (±0.03)	x
	Weak layer thickness (m)	0.015 (±0.002)	x	x	x	0.11 (±0.01)		0.21 (±0.02)
	Substratum thickness (m)	0.275 (±0.03)				0.12 (±0.01)		1.27 (±0.05)
	Esl (MPa)	5 (±0.5)	x	x	x	9 (±3)	x	5 (±6)	x
	Ewl (MPa)	0.41 (±0.09)		x	x	0.41	x	0.41	x
	Esub (MPa)	13(±0.09)			x	90 (±60)		50 (±40)
	Poisson's ratio	0.25 (±0.05)		x	x	0.25 (±0.05)	x	0.25 (±0.05)	x
	Density slab (kg m -3 )	159 (±8)	x	x	x	181 (±9)	x	157 (±44)	x
	Density weak layer (kg m -3 )	138 (7)	x		x	-		-
	Density substratum (kg m -3 )	325 (±58)			x	296 (±44)		265 (±46)
	Collapse height (mm)	9.3 (±0.9)		x		2 (±0.6)	x	4(±2)
	2.2 Numerical models
	2.2.1 Discrete Element Method

Bobillier and others (2020)

presented a three-dimensional discrete element method (DEM) model to reproduce failure in a layered snow sample. DEM, first introduced by Cundall and Strack

Table 2 : Theoretical wave speeds for the PST, whumpf and avalanche. The longitudinal wave speed 𝒄 𝒍 , the shear wave speed 𝒄 𝒔 and Rayleigh wave speed 𝒄 𝒓 of the slab was computed from snow properties. The three columns on the right show crack propagation speed estimates. To calculate 𝒄 𝒔𝒄 and 𝒄 𝒇𝒘 , the formulations by McClung (2005) and Heierli (2005), respectively, were used. In the last column the experimental crack speeds 𝒄 𝒆𝒙𝒑 are listed for comparison.

 2

			Elastic wave speeds		Crack speed estimates
		𝑐 𝑙 (m s -1 )	𝑐 𝑠 (m s -1 )	𝑐 𝑟 (m s -1 )	𝑐 𝑠𝑐 (m s -1 )	𝑐 𝑓𝑤 (m s -1 ) 𝑐 𝑒𝑥𝑝 (m s -1 )
	PST	177 ±10	112 ±7	103 ±6	(79 -101) ±6	35.8 ±1.4	39 ±13
	Whumpf	224 ±35	105 ±6	97 ±6	(74 -95) ±6	54 ±6	49 ±5
	Avalanche	170 ±90	152 ±34	140 ±31	(106 -137) ±30	39 ±11	36 ±6

Acknowledgements. We thank Laurent Crettenand for capturing the video of the avalanche. We are also grateful to Augustin Rion and Armand Salamin who recorded the corresponding snow profile. Further, we thank the reviewers, Pascal Hagenmuller and Karl Birkeland for their positive and constructive comments that helped to improve the manuscript. BB was supported by a grant from the Swiss National Science Foundation (200021_169424).

References:

Crack speed estimates from the avalanche video may be subject to systematic errors. To estimate crack speed, we tracked slab fractures that progressively became visible on the video images. With this method, we thus assumed that slab fractures are a good proxy for the location of the crack tip in the weak layer. However, there are three processes that can result in a lag between the appearance of the slab fracture and the actual location of the crack tip, and thus consistently lower crack speed estimates. First, at the trigger point, the explosive likely destroys the weak layer over an area that extends several meters beyond the explosion. Thus, crack propagation does not start at x = 0 m but rather at x = x0 ( a). Second, a slab fracture may only start to form when the crack in the weak layer has already propagated a distance x1 further. Third, there is a delay between the formation of the slab fracture and the first time it becomes visible on the images (x2). Since the contributions x0 and x1 are unknown, we assumed these to be constant and combined them 𝑥 ′ = 𝑥 0 + 𝑥 1 (Fig. 1). To estimate x2, we assume that a slab fracture becomes visible when it reaches a size of one full pixel on the image 𝑥 𝑣𝑖𝑠 = 1 𝑝𝑖𝑥𝑒𝑙. The time 𝑡 𝑣𝑖𝑠 needed for a slab fracture to open this distance can be estimated as:

with 𝑔 the gravitational acceleration, 𝜃 = 39° the mean slope angle and 𝜇 = 0.57 the frictional resistance [START_REF] Van Herwijnen | Measurement of crack-face friction in collapsed weak snow layers[END_REF] after weak layer cracking.

Once the slab fracture opening time 𝑡 𝑣𝑖𝑠 is estimated, the error 𝑥 2 is estimated as 𝑥 2 = 𝑐 𝑡 𝑣𝑖𝑠 , where c is the crack speed in the weak layer. Finally, we assume constant crack speed in the weak layer and interpreted our measured data pairs 𝑡 𝑚 , 𝑥 𝑚 with a simple model accounting for the uncertainties 𝑥 2 and 𝑥 ′ :

where 𝑡 𝑚 and 𝑥 𝑚 are the initial estimates for crack propagation time and distance leading to the crack speed measures (green squares in Fig. 1). First, we neglected the contributions of 𝑥 ′ , and just optimized (in a least squares manner) crack speed c in Equation C1 with our data pairs 𝑡 𝑚 , 𝑥 𝑚 (Fig. 1b, orange line). This results in a weak layer crack speed of c = 43 ± 8 m s -1 . 705 Second, we also accounted for the error 𝑥 ′ (Fig. 1b, blue line). This additional degree of freedom in the model had little influence on the crack speed (c = 45 ± 5 m s -1 ) and 𝑥 ′ was estimated to be 𝑥 ′ ≈ 15 m. This simple model provides a plausible explanation for the observed initial increase in crack speed. Still, the assumption of a terminal crack speed cannot be validated. For larger distances, this assumption is supported by our data, as crack speed estimates were relatively constant, and this is the range where the relative influence introduced by the error contributions 710 becomes negligible (e.g. relative deviation against terminal speed c(x > 200m) ≤ 5 %).