Influence of slab depth spatial variability on skier-triggering probability and avalanche size

Abstract

Spatial variability of snowpack properties adds uncertainty in the evaluation of avalanche hazard. We propose a combined mechanical–statistical approach to study how spatial variation of slab depth affects the skier-triggering probability and possible release size. First, we generate multiple slab depth maps on a plane fictional slope based on Gaussian Random Fields (GRF) for a specific set of mean, variance and correlation length. For each GRF, we derive analytically the Skier Propagation Index (SPI). We then simulate multiple skier tracks and computed the probability based on the number of skier hits where SPI is below 1. Finally, we use a depth-averaged material point method to evaluate the possible avalanche size for given slab depth variations. The results of this analysis show that large correlation lengths and small variances lead to a lower probability of skier-triggering as it reduces the size and the number of areas with low slab depth. Then, we show the effect of skiing style and skier group size on skier-triggering probability. Spatial variability also affects the possible avalanche size by adding stress fluctuation causing early or late tensile failure. Finally, we demonstrate with our models the well-known relationship between the probability and the size in avalanche forecasting.

Introduction

Avalanche hazard

Snow avalanches represent a natural hazard for infrastructure and backcountry recreationists across mountainous areas all across the world (Stethem and others, 2003). Snow avalanches can be divided into different types from wet or dry avalanches and also loose-snow or slab avalanches (Schweizer and others, 2003). However, dry-snow slab avalanches are the most destructive type and also very difficult to predict (Techel and others, 2016). Such an avalanche releases from the failure of a porous weak layer buried below a cohesive snow slab. The initiation of the failure can be induced by a skier, new snow or explosives. If the size of the failed zone in the weak layer – or crack – exceeds a critical size, the crack may self-propagate across the slope possibly leading to the release and sliding of the snow slab (Schweizer and others, 2016). Statham and others (2018) proposed a conceptual model to better predict the avalanche hazard for practitioners and backcountry recreationists. This conceptual model consists of two main components: the likelihood of an avalanche and the avalanche release size. The combination of these two components gives the avalanche hazard in North America (Statham and others, 2018). The terminology is slightly different in Europe which is based on three components: the probability of avalanche release, the frequency of these triggering spots and the avalanche size (Techel and others, 2020). Practitioners and forecasters use mostly snow stability tests to gain information on the probability of triggering an avalanche, both by a skier or naturally. Repeated tests at numerous locations help them to gain information on the spatial distribution of the instability (Schweizer and others, 2008). However, the sparse and punctual nature of available observations on snowpack properties makes the forecasting of dry snow slab avalanches difficult (Hägeli and McClung, 2004). The second main component of the avalanche hazard is the avalanche release size. Forecasters try to estimate the possible avalanche size by estimating the slab depth and the slab-weak layer propagation propensity. Studies have shown that these properties are highly variable in space (Schweizer and others, 2008), which makes it even more difficult to predict. The snow spatial variability at different scales also adds complexity to this difficult task by adding uncertainty on whether the properties measured on the field are representatives or not of the slab and weak layer system (Schweizer and others, 2008). At a smaller scale, the decision-making in avalanche terrain such as the task of up-hill and downhill route finding is very complex partly due to the spatial variability of snow.

Spatial variability of snow

The subject of the spatial variation of snow mechanical properties is nothing new in the science community. Conway and Abrahamson (1984) has made several measurements of a weak layer shear strength along an avalanche fracture line. They highlight these so-called deficit zones where the weak layer shear strength was significantly lower than surrounding areas along the fracture line. The idea of weak zones or weak spots initiated numerous studies over the last two decades (Schweizer and others, 2008). New instruments like the high-resolution snow penetrometer called SnowMicroPen (SMP) can measure with great accuracy the snow mechanical properties of both the slab and the weak layer (Johnson and Schneebeli, 1999; Löwe and van Herwijnen, 2012; Proksch and others, 2015; Reuter and others, 2019). The SMP enabled fast sampling of the snow mechanical properties over several locations on an avalanche-prone slope (Landry and others, 2004; Kronholm, 2004; Feick and others, 2007; Lutz and others, 2007; Bellaire and Schweizer, 2011; Lutz and Birkeland, 2011; Reuter and others, 2016). These studies focused mainly on the weak layer properties, thickness and strength, and they used the variogram analysis to estimate the spatial pattern of the weak layer mechanical properties. They reported correlation distances ranging from 1 to 15 m with a spatial sampling extent of around 20 and 40 m. The correlation distances were nearly half of the extent and this could bias the estimation of the correlation length (Skøien and Blöschl, 2006; Kronholm and Birkeland, 2007). Reuter and others (2016) used the same sampling density but with a spatial extent of nearly 500 m, and reported correlation distances from 5 to 25 m with one exception of 68 m.

These studies measured and explained the spatial variability of snow mechanical properties, but did not explain the effect of this variability on slope stability. Kronholm and Schweizer (2003) proposed a conceptual model to explain the effect of the spatial variation of the stability on the overall slope stability. Short-range variation could have a stabilizing effect on the snowpack and the long-range variation could have a so-called ‘knock-down effect’ on the slope stability, but further investigation through mechanical models was needed to test this conceptual model. Several studies simulated artificial spatial patterns of the weak layer into mechanical models to explain the effect of the spatial variability of the weak layer on the overall slope stability (Fyffe and Zaiser, 2004; Kronholm and others, 2004; Schweizer and others, 2008; Gaume and others, 2013, 2014). First, some studies used cellular automata models and showed the effect of the variance shear strength on the slope stability (Fyffe and Zaiser, 2004; Faillettaz and others, 2004; Kronholm and Birkeland, 2005). High shear strength variances create more deficit zones and cause easier overall failure of the slope even with strong zones. However, these models only account for the state of the neighbor cells and large-scale elastic redistribution could not be taken into account, therefore the link with the correlation length could not be explored. Finite-element method (FEM) can handle large-scale elastic redistribution and was used in several studies to explain the influence of the correlation length over the slope stability. Gaume and others (2013) explored how the stress redistribution induced by the elasticity of the slab could smooth the heterogeneity of the weak layer. They showed, in particular, that if the correlation length is smaller than the characteristic elastic length of the system, it behaves as in a homogeneous case. When the correlation length is larger than this elastic length, the smoothing does not take place and the system is more likely to fail even for large slab depth which illustrated the so-called knock-down effect (Gaume and others, 2013, 2014). Gaume and others (2015) use the same method to estimate the propensity for tensile failure in the slab in order to relate to the avalanche release size. Weak layer heterogeneity increases slab tensile failure propensity for soft and shallow slabs, thus potentially smaller avalanches. On the other hard and deep slabs were hardly influenced by weak layer heterogeneity which led to wide-spread crack propagation.

These FEM studies focused exclusively on the spatial variation of the weak layer cohesion and its influence on natural avalanche release. It has been demonstrated numerous times, through survey and modeling, that snow depth is highly variable in mountainous areas (e.g. Winstral and others, 2002; Deems and others, 2006; Mott and others, 2011; Schirmer and others, 2011; Grünewald and others, 2013; Hubbard and others, 2018), especially in avalanche start zones (Miller and others, 2022). The slab depth spatial variability should be related to the spatial variation of snow depth. Therefore, the spatial variation of the slab depth should affect the skier-triggering probability over an entire slope. For homogeneous cases, the influence of slab depth on skier-triggering was analyzed based on stress over strength approaches (Föhn, 1987; Habermann and others, 2008; Monti and others, 2016) but also to address crack propagation propensity (Heierli and others, 2008; Gaume and others, 2017). To the best of our knowledge, there is no study on the effect of slab depth spatial variation on the probability of skier-triggering and crack propagation. Slab depth variability should also affect the avalanche possible size. Slab properties, mainly slab depth and density, are one of the main drivers for dynamic crack propagation (McClung, 1981; Heierli and others, 2008; Gaume and others, 2017). The spatial variation of the slab-weak layer system should affect the stress in the slab and might promote crack arrest caused by slab tensile failure. There is a need for further investigation on that matter in order to provide estimates of the avalanche release size which is crucial information for avalanche forecasting and risk management.

Several tools have been used to understand dynamic crack propagation. Previous studies on the influence of spatial variability of snow on the release and the avalanche size used mainly mesh-based approaches which can handle large deformations and fracture propagation at the (significant) cost of re-meshing and mesh refinement techniques. Some researchers used the discrete element method to model snow but are limited to small computational domains from the microstructure scale (Hagenmuller and others, 2015; Mede and others, 2020) to a few meters scale (Gaume and others, 2017; Bobillier and others, 2021). More recently, the material point method (MPM), a hybrid Eulerian–Lagrangian technique showed promise in simulating solid–fluid transitions and crack propagation in geomaterials (Sulsky and others, 1994). Gaume and others (2019) propose an elastoplastic MPM with a new snow constitutive model to simulate the mechanical behavior of snow slabs and weak snow layers. With this model, they can simulate in a unified way, numerous important mechanical processes in a snow slab avalanche, from the weak layer failure initiation, dynamic crack propagation, slab tensile fracture, and, finally the release and the flow of the slab down the slope. Trottet and others (2022) used this method and full-scale measurements to reveal a transition between anticrack propagation (closing crack/mode – I) and so-called ‘supershear’ propagation in which the fracture occurs in shear (mode II) at a speed larger than the shear wave speed (of the slab). They report that this transition requires a slope angle larger than the snow friction angle (≈ $27^\circ$) and a propagation distance larger than typically 3–5 m. This suggests that a pure shear interface model could be sufficient to simulate large avalanche releases on an inclined slope. This finding, together with the large computational cost of 3D MPM simulations motivated us to use the depth-averaged MPM (DAMPM) recently proposed and validated by Guillet and others (2023). DAMPM enables a fast computation time over large domains, allowing us to produce a vast number of simulations to better understand the effect of the slab depth variation on the avalanche release size.

Understanding how the spatial variation of slab depth influences the probability of skier-triggering and the avalanche release size is of great importance. This could lead to an improvement in snow avalanche forecasting and decision-making in avalanche terrain for backcountry recreationists. To better understand this issue, we propose a combined mechanical–statistical approach to study how spatial variation of slab depth affects the skier-triggering probability and possible avalanche release size. First, a sensitivity analysis is made for the skier-triggering probability in relation to the spatial variation of the slab depth. Secondly, we use our approach to test the effect of skiing style and downhill strategy on skier-triggering probability. Lastly, DAMPM simulations allowed us to investigate the relation between slab depth variations and the avalanche release size, by analyzing the propagation distance leading to the first tensile fracture in the slab.

Methods

Spatial variability of slab depth

Gaussian random fields (GRF) will be used to generate artificial 2D surfaces as input into our numerical mechanical models. Each point of the surface in space (x,  y) can be defined as a random variable, so the collection of random variables is called a random field defined by:

(1)$$D( x,\; y) \in \overline{D} \cdot C$$

where D(x,  y) is the slab depth random field, $\overline {D}$ is the mean slab depth and C is the Gaussian covariance model defined by the formulation in gstools v1.3 documentation (Müller and others, 2022):

(2)$$C( d) = S^2_D \left[1- \exp \left({d\over \epsilon}\right)^2 \right]$$

where d is the distance between the observations, $S^2_D$ is the slab depth variance and $\epsilon$ is the correlation length of the slab depth variation. We used a Gaussian covariance model without a nugget. For each GRF, we can specify the mean and the covariance function, the last one is defined by a variance and a correlation length. A sensitivity analysis is presented on these three parameters of the slab depth spatial variation for both the skier-triggering probability and avalanche release size. We test the mean slab depth from 0.5 m to 1 m with a 0.1 m increment. Several studies have shown that snow mechanical properties can be approximated by a Spherical and Gaussian covariance function with correlation lengths ranging mostly from 5 to 20 m (Kronholm and Schweizer, 2003; Landry and others, 2004; Feick and others, 2007; Lutz and others, 2007; Bellaire and Schweizer, 2011; Lutz and Birkeland, 2011; Reuter and others, 2016) and sometimes 50 m (Reuter and others, 2016). Based on these studies, we chose to test correlation length from 5 m to 40 m with an increment of 5 m. The idea is to tend toward a very smooth random field close to a homogeneous case. We took the same approach with the variance, from 0 to 0.025 m² which represents approximately 0.3 m from the mean slab depth to maximum and minimum values. The selected variance values were based on field measurements done by Meloche and others (2023). Gaussian random fields were generated using the gstools v1.3 package in Python (Müller and others, 2022).

Skier-triggering probability

One of the specific objectives of this study is to evaluate the skier-triggering probability as a function of slab depth variations. The snow properties were selected to represent a meta-stable slab and weak-layer system where the natural release is not possible but only the skier-triggering. Different snow mechanical properties can be related to each other (Jamieson and Johnston, 1990; Scapozza and others, 2004; Sigrist, 2006). In general, it has been demonstrated that the increase in slab depth will also increase the snow density, elastic modulus, and shear strength of the weak layer. We computed realistic snow mechanical values based on empirical power-law functions obtained based on literature data. First, the snow density ρ is related to the mean slab depth $\bar D$ according to Eqn (6) in McClung (2009):

(3)$$\rho = 100 + 135 \bar D^{0.4}$$

Then, we computed the elastic modulus E based on the snow density ρ according to Sigrist (2006):

(4)$$E = 9.68 \times 10^8 \left({\rho\over \rho_{ice}}\right)^{2.94}$$

where ρ _ice = 917 kg m⁻³. Note that, contrary to the slab depth, the density and elastic modulus of the slab does not spatially vary in our simulations but are related to the mean slab depth. Finally, the shear strength of the weak layer τ _p was computed according to the following power–law relationship from Bažant and others (2003) adapted by Gaume and others (2014):

(5)$$\tau_p = c + 1370D^{1.3}$$

where c is the cohesion set to 300 Pa. This allows us to account for the local friction effect where the slab depth is locally increasing.

Stability metrics

We choose to assess the skier stability based on the skier propagation index (SPI) proposed by Gaume and Reuter (2017). The index consists of the ratio between two lengths, namely the critical crack length a _c and the skier crack length l _sk (Gaume and Reuter, 2017). The skier crack length is computed by solving the equation: τ + Δτ = τ _p, where $\tau = \rho g D \sin \psi$ is the shear stress due to the slab weight on an inclined slope ψ and g is the gravitational acceleration constant. The additional stress due to the skier standing on the snow is defined by (Föhn, 1987; Monti and others, 2016):

(6)$$\Delta \tau = {2 R \cos \alpha \sin^2 \alpha \sin( \alpha + \psi) \over \pi D}$$

where R is the skier load set to 780 N and ψ is the slope angle set to $35^\circ$. We find the two roots of the equation given two angles α ₁ and α ₂, where total shear stress is equal to the shear strength of the weak layer. We used the two angles α ₁ and α ₂ to find the skier crack l _sk with the following equation:

(7)$$l_{sk} = D \left[{1\over \tan \alpha_1} - {1\over \tan \alpha_2}\right]$$

The critical length a _c is computed using the formulation from Gaume and others (2017):

(8)$$a_c = \Lambda \left[{- \tau + \sqrt{ \tau^2 + 2\sigma ( \tau_p - \tau) }\over \sigma}\right]$$

where $\sigma = \rho g D \cos \psi$ and Λ is a characteristic length of the system defined by:

(9)$$\Lambda = \sqrt{{E^{'}D D_{wl}\over G_{wl}}}$$

with E ^′ = E/(1 − ν ²). The weak layer thickness D _wl and the shear modulus G _wl were set to 0.04 m and 0.2 MPa respectively, ν is the Poisson ratio set to 0.3. The skier propagation index SPI is thus defined as:

(10)$$SPI = {a_c\over l_{sk}}$$

These formulas allow us to generate SPI maps from our corresponding GRF slab depth and snow properties maps over a 200 m by 100 m fictional slope (Fig. 1). For each GRF realization, we simulated sinusoidal ski tracks all over the fictional slope representing a ‘modern freeride’ skiing trajectory defined by a down-slope turn radius of 10 m and cross-slope amplitude of 5 m (Fig. 1). The spacing between the skier was held constant at 5 m with a total of 40 skiers. We recorded a hit if a skier track passed through a zone with SPI below 1. We computed the probability of skier-triggering with the number of skier tracks who recorded a hit compared to all skier tracks on the slope. A convergence of analysis for the total number of skiers on the slope is presented in the appendix (Fig. 11).

Figure 1. Example of one realization of a Gaussian Random Field GRF for slab depth D, slab density ρ, skier crack length l _sk, critical crack length a _c, skier propagation index SPI and skier tracks with SPI. The slope is reduced to 50 m by 100 m only for visualization purposes.

A secondary objective of this study was to assess the influence of skiing style or skiing trajectory on the skier-triggering probability. Different trajectories were tested on the basis of two ‘extreme’ trajectories. The first one is to mimic a pure linear trajectory which is defined by the projected weak spot length on the x-axis (cross-slope) compared to the total cross-slope length. This ratio was obtained following these two steps: (1) the presence of weak spots in the up-slope direction is checked for every point following a transect in the cross-slope direction, and (2) a length is obtained with the sum of every presence of weak spot along this transect and then compared to the cross-slope length. The second extreme trajectory mimics a skier who will span the entire slope in a cross-slope direction which would yield a skier-triggering probability of one.

Possible avalanche size

To link the spatial variability of slab depth to the avalanche release size, a method is needed to compute the dynamic crack propagation in a weak snow layer. The finding from Trottet and others (2022) supports the physical assumptions needed for a depth-averaged method that integrates the momentum and mass conservation equation in the z direction. Guillet and others (2023) presents a detailed view of this method and the integration of the governing equation. We will present the key assumptions but please refer to Guillet and others (2023) for an in-depth view of the method. The first key assumption is the classic shallow water assumption where the vertical length is shorter than the horizontal length, which translates in an avalanche context as ${h_0\over L} \ll 1$ where h ₀ is the standard height of the slab and L the characteristic length of the avalanche. The material is assumed incompressible meaning that the density ρ does not depend on position either time and place, and the flow surface is stress-free at the top of the slab where σ|_h=z = 0. Note here, that the method could be easily adapted in compressible form if the ρ needs to be as a function of the particle height and the position. The velocities in the x direction are similar to the ones in the y direction but the velocities in the z direction are small and negligible. The other particularity of the depth-averaged method is the column-shaped particle of height h = D where the integration point of σ _zz is h/2 (Fig. 2), which gives us this depth-averaged equation of mass conservation:

(11)$${\partial h\over \partial t} + {\partial ( h \overline{v}_x) \over \partial x} + {\partial( h\overline{v}_y) \over \partial y} = 0$$

where $\overline {v}_x$ and $\overline {v}_y$ are the depth-averaged velocity field at the integration point of the particle. The last assumption is a plane stress assumption with σ _zy = σ _zx = 0. Basal forces are applied at the interface between the slab and the weak layer defined by:

(12)$$\tau_{xz}\, \colon = \sigma_{xz}\vert _{z = 0},\; \tau_{yz}\, \colon = \sigma_{yz}\vert _{z = 0}$$

which give also σ _zz|_z=0 = ρgh at the bottom of the slab column and $\overline {\sigma }_{zz} = {1\over 2}h\rho g$ at the integration point. This allows us to have the depth-averaged non-conservative form of momentum conservation:

(13)$$\rho h {\overline{d}\overline{v}_x\over \overline{d}t} = {\partial( h \overline{\sigma}_xx) \over \partial x} + {\partial( h \overline{\sigma}_xy) \over \partial y} - \tau_{xz} + g \rho h .$$

(14)$$\rho h {\overline{d}\overline{v}_y\over \overline{d}t} = {\partial( h \overline{\sigma}_yy) \over \partial y} + {\partial( h \overline{\sigma}_xy) \over \partial x} - \tau_{yz} + g \rho h .$$

Figure 2. DAMPM model and simulation parameters. (a) The geometry of the simulated propagation saw test PST with a sinusoidal slab depth variation and shear strength associated. The slab depth variation parameters from the sinus function are the mean slab depth $\overline {D}$, the standard deviation S _D, and the correlation length $\epsilon$. (b) Slab elastoplasticity model following a Modified Cam-Clay yield surface. (c) Weak-layer model as quasi-brittle interface.

The slab is represented with an elastoplastic model following a Modified Cam-Clay yield surface γ in the p − q space (Gaume and others, 2018):

(15)$$\gamma( p,\; q) = q^2( 1 + 2\beta) + M^2( p + \beta p_0) ( p-p_0) $$

where M is the cohesionless critical state line, β is the cohesion parameter that quantifies the ratio between the tensile and compressive resistance and p ₀ is the consolidation pressure that affects the size of the yield surface. With M and β being constant respectively at 0.3 and 1.2, we related the slab tensile strength σ _t to p ₀:

(16)$$p_0 = -( 1-\beta) + \sqrt{\left({( \beta + 1) p\over 2\beta}\right)^2 + \left({q\over M}\right)^2 \left({( 1 + 2\beta\over \beta}\right)}$$

where $p = {\sigma _t\over \sqrt {3}}$ and $q = {\sigma _t\over 2}$ considering the plane stress hypothesis. The weak layer is an external force represented by a quasi-brittle interface with a softening behavior when the displacement u reaches the critical displacement u _c defined by the shear stress reaches τ _p (Fig. 2). Then we have the softening until the displacement reaches the residual displacement u _r defined by the residual friction τ _r with the snow friction angle Φ = 27°. The softening displacement is set to δ = 2u _c, as suggested in Gaume and others (2013). The snow parameters are the same as in the skier method regarding snow density and elastic modulus, as they are set according to the mean slab depth values of the simulated slope via a power–law fit. The shear strength is set according to the local slab depth values to represent the friction for thicker slabs. The slope angle is also the same as the previous method set at $35^\circ$. We used the empirical power–law fit to link the slab tensile strength with the mean slab density according to Sigrist (2006):

(17)$$\sigma_t = 2.4 \times 10^5 \left({\rho\over \rho_{ice}}\right)^{2.44}$$

The depth-averaged material point method is more computer-intensive and time-consuming compared to the analytical method. A GRF simulation of the slab depth required at least 50 simulations per set of GRF parameters to obtain a statistical distribution and get representative results. We changed our approach and simulated a sinusoidal slab depth variation to reduce the number of simulations to one and have a deterministic view. The variance parameter is changed to standard deviation S _D representing the amplitude of a sinus function (Fig. 2). We simulate very long propagation saw tests (PST) of 75 m in the up-slope direction and 0.30 m wide (Gauthier and Jamieson, 2008). For each simulation, the cohesion is numerically removed from the bottom of the slope to obtain a critical crack length that will start the dynamic crack propagation across the slope. Finally, a tensile fracture is needed in the slab to be released from the slope and create a slab avalanche. This type of fracture in the slab occurs when the tensile stress in the down-slope direction (σ _xx) reaches the tensile strength of the slab σ _t. The distance to the first tensile fracture in the slab is noted for each simulation, we called this distance the tensile length L _t.

Results

Sensitivity analysis of the skier-triggering probability

The probability of skier-triggering decreases with the increase of the mean slab depth as expected (Fig. 3). The homogeneous cases with zero variance show stable snowpack for skiers triggering from 0.6 m to 1 m slab depth. We observed two distinct regimes, one for a slab depth of 0.5 m where every skier is triggering for both homogeneous or heterogeneous cases. The second regime is defined with a 1 m mean slab depth where every skier doesn't trigger except for the larger variance at 0.025 m². There is an intermediate regime for the other values of mean slab depth which is affected by both the variance and the correlation length. The probability of skier-triggering increases with the slab depth variance. The increase in the slab depth variance allows some areas of the fictional slope to have shallow slab depth and create weak zones where the SPI is below 1. The increase of the correlation length is decreasing the skier-triggering probability. This decrease is caused by the reduction of the number of weak spots with SPI below 1 across slopes. A small correlation length creates more weak spots where SPI is below 1 and a large correlation length creates fewer but bigger weak spots with SPI below 1. Instead of the mean presented above, we looked into the distribution of the skier-triggering probability containing 100 realizations. Figure 4 shows different probability density functions for each specific set of GRF parameters. With a mean probability below 0.3, the distribution is skewed towards 0 like a log-normal distribution. For cases where the mean skier-triggering probability range from 0.3 to 0.6, the distribution is relatively flat and centered, like a normal distribution. Even if the mean value is around 0.4 or 0.5, the variance of these distributions is quite large from 0.2 to 0.8 skier triggering probability. This means that for the same GRF parameters, some realizations were very unstable with more than half of the skiers triggered an avalanche and some realizations with very few skiers triggered. These cases in the intermediate regime described above from 0.6 to 0.9 m mean slab depth, the spatial variability of slab depth adds uncertainty to the probability to trigger an avalanche by a skier. With a mean probability above 0.8, the distribution is normally distributed around the mean values.

Figure 3. Sensitivity analysis for the mean skier-triggering probability with regards to the mean, spatial variance and correlation length of the slab depth. The mean probability is computed from a distribution of 100 realizations for each specific set of mean slab depth, variance and correlation length.

Figure 4. Probability density function of the skier-triggering probability for mid to low mean skier-triggering probability. All the distributions presented are from a GRF using a 0.7 m mean slab depth, 0.0025–0.0075 m² slab depth variance, and 10-20-30 m correlation length.

Influence of the skiing style

The previous analysis was conducted using a constant skiing style defined by a 10 m down-slope turn radius (R _down) with a 5 m cross-slope amplitude (A _cross) (Fig. 1). However, these two parameters should influence the trajectory of each skier and therefore, should also affect the probability of the skier hitting a weak spot. We made a sensitivity analysis of these two parameters. We test multiple values to get a down-slope linear trajectory and the complete opposite with a trajectory that will traverse the entire slope in the cross-slope direction. These two extreme trajectories should mimic a freeride down-hill skiing trajectory (linear down-slope) compared to an up-hill skinning trajectory (cross-slope). Figure 5 shows that for a large A _cross/R _down ratio, the mean skier probability is increasing towards a probability of one with positive hits for every skier, meaning that a trajectory with a large cross-slope amplitude is more likely to encounter weak spot and trigger an avalanche. For a small A _cross/R _down ratio, the mean skier probability is decreasing towards a value determined by a ratio that we called the linear weak spot cross-slope ratio (${\sum {\hat {L}_x}}/{L_x}$). This ratio indicates the minimum probability of a purely linear trajectory specific to each realization per set of GRF parameters. Note that the 0.505 found in (Fig. 5) represents the mean of the linear weak spot cross-slope ratio for 100 realizations for a mean slab depth of 0.7 m, variance of 0.0075 m², and 20 m correlation length.

Figure 5. Mean skier-triggering probability from 100 realizations for different skiing style ratios A _cross/R _down. A _cross represents the cross-slope amplitude and R _down represents the down-slope turn radius. A small skiing style ratio represents a linear down-slope trajectory and a large skiing style ratio represents a cross-slope trajectory. The probabilities are constrained by two values: the first is a probability of 1 where all skiers have triggered, and the second value is a linear weak spot cross-slope ratio between the sum of weak spot length in the cross-slope direction compared to the total cross-slope length. The dashed line is set at 0.505 which is the mean of the linear weak spot cross-slope ratio for 100 realizations for this set of GRF parameters: mean slab depth of 0.7 m, variance of 0.0075 m², and 20 m correlation length. This line should move with regard to the GRF parameters. The inlets represent schematic skiing style based on the skiing style ratios.

Influence of group size and terrain choice

We investigated the influence of the group size on the skier-triggering probability. Here, we want to test the hypothesis that it is safer to ski near a preexisting ski track that is considered safe, compared to a completely random approach. We randomly selected a position at the top of the slope to start the first skier, and if the first skier track didn't record a hit on a weak spot, we keep adding a skier 5 m apart until a trigger was recorded or the slope was entirely skied. We repeated this operation 50 times on one realization and repeated it for 100 realizations for a 0.7 m mean slab depth and 0.005 m² variance, with four different correlation lengths (5, 10, 20, 30 m). Figure 6 is showing that for a correlation of 5 and 10 m, the random and structured approaches have a similar experimental cumulative density function (ECDF) but the structured approach is slightly shifted. The ECDF for the 5 and 10 m correlation length had a median of one or two additional skiers which can be translated into a slightly lower skier-triggering probability compared to the 20 and 30 m ECDF. For the simulation with a correlation length of 20 and 30 m, the ECDF of the structured approach is completely shifted from the random approach, towards more additional skiers before triggering. The ECDF is starting at two additional skiers for the 20 m and the 30 m distribution, which means that a minimum of two additional skiers were needed before recording a trigger compared to the random approach, from all the 5000 simulations represented in one ECDF. The difference is more significant for the median of the ECDF, the random approach as between 3 to 4 additional skiers before the trigger compared to the structured approach with 7 skiers (20 m) and 10 skiers (30 m).

Figure 6. Experimental cumulative density functions (ECDF) of the number of additional skier to trigger from mid to low mean probability. All the distributions presented are from a GRF using a 0.7 m mean slab depth, 0.005 m² slab depth variance, and 5-10-20-30 m correlation lengths. The inlets show examples of the corresponding SPI map specific to their GRF parameters.

Tensile length and avalanche release size

Figure 7 shows different simulations of PST for three different correlation lengths. During a dynamic crack propagation for a homogeneous case, the tensile stress σ _xx increased linearly as the unsupported slab grew from the damaged weak layer. Then the tensile stress σ _xx reached the tensile strength σ _t, and a tensile fracture in the slab is observed. The increase is linear because the shear stress is constant across the slope in a homogeneous case, as the crack tip moves across the slope at 1.6 C _s. However, the tensile stress pattern changed for different cases of slab depth variations. If the correlation length is large, around 25 m or more, the tensile length converges towards the homogeneous case because the tensile stress is increasing almost constantly. Figure 7a shows a simulation for a mean slab depth of 0.7 m, a standard deviation of 0.25 m, and a correlation length of 30 m where the tensile stress σ _xx is constantly increasing in the lower part of the slope (distance up-slope < 20 m). Then, the tensile stress increases non-linearly to reach the tensile strength at ≈ 34 m. Tensile stress σ _xx increases rapidly as shear strength τ _p and shear stress τ _xz reduces around 25 m. This smooth spatial variability with a long correlation length results in a tensile length similar to the homogeneous case. Figure 7b shows a simulation for a mean slab depth of 0.7 m, a standard deviation of 0.25 m, and a correlation length of 15 m where the tensile stress fluctuates as the crack tip moves across the slope. The same nonlinear increase is observed around 15 m up-slope in the tensile stress, but then it decreases as the slab depth, the shear stress, and strength increase again around 25 m up-slope (Fig. 7b). Finally, the tensile stress increases a second time as the slab depth, shear stress and strength start to decrease around 40 m up-slope (Fig. 7b). In this particular case, the spatial variability of the slab depth and the underlying weak layer strength cause a fluctuation in the tensile stress of the slab σ _xx, resulting in a longer tensile length. Figure 7c shows the same fluctuation in tensile stress but for a 10 m correlation length. The fluctuation in tensile stress is more significant and closer to the bottom of the PST, leading to a shorter tensile length. From a static point of view, tensile stress builds as the length of the ‘unsupported’ slab increases due to weak layer damage and crack propagation. The tensile stress is equal to the load from the slab in the down-slope minus the friction of the slab weight. This relation is linear assuming a constant slab depth. However, our system had a variation in slab depth that will also cause a variation in the friction effect that is more pronounced where the slab is locally thicker and explains why the tensile stress locally reduces where the slab is thicker (Fig. 7). Dynamic effects could also explain some fluctuations in the tensile stress. In fact, variations in shear stress and weak layer strength lead to variations in crack speed around 1.6 C _s, which can also cause fluctuations in tensile stress.

Figure 7. Three different regimes of tensile fracture for a 0.7 m slab depth $\overline {D}$ and 0.25 m S _D. These three simulations show the last frame saved just before the tensile fracture occurs when σ _xx = σ _t. The shear stress τ _xz and the weak layer shear strength τ _p and the cohesion c are also represented. The crack tip (dotted line) is located just behind the peak of τ _xz at the loss of cohesion c. The distance between the crack tip and the τ _xz peak is due to the softening δ. The bottom plot shows the corresponding crack speed $\dot {a}$ which is normalized over the shear wave speed $C_s = \sqrt {{G}/{\rho }}$, and the slab depth D in m. (a) PST simulation with a $\epsilon$ 30 m just before a tensile fracture occurs. (b) PST simulation with a $\epsilon$ of 15 m just before a tensile fracture occurs. (c) PST simulation with a $\epsilon$ of 10 m just before a tensile fracture occurs.

We present a sensitivity analysis of the propagation distance that results in the first tensile fracture (L _t) from our numerical PST experiment. Figure 8 shows all the results from our simulations with different mean slab depth values $\overline {D}$, slab depth standard deviation S _D (sinus function) and correlation length $\epsilon$. The homogeneous tensile length is obtained here based on a simulation with homogeneous slab depth (i.e. standard deviation of zero). This homogeneous tensile length increased with the slab depth even if the theoretical quasistatic tensile length is in principle not linked to the slab depth ($L_t = {\sigma _t\over \rho g \sin \psi }$). However, in our study, the slab density and tensile strength are related to the mean slab depth Eqns (3) and (4), which explains the reported increase in the tensile length obtained for our homogeneous cases (Fig. 8). A shorter tensile length is obtained as the standard deviation increases. A longer tensile length is obtained as the correlation length increases around 25 m, then, it converges towards the same values as the homogeneous case. We observed a tensile length approximately 20 m shorter than the homogeneous case when the standard deviation is 0.2 m and a correlation length of 10 m or less. The tensile length is only longer than the homogeneous case for a correlation length of 15 m and 20 m for higher slab depth, but also for a standard deviation of approximately 0.1 m and higher (Fig. 8). This regime with larger tensile lengths is more pronounced for thicker slab depth. This particular regime with longer tensile lengths, which is associated with a correlation length of 15 or 20 m, is caused by the use of the sinus function. This function caused an increase in the weak layer strength around 30–40 m, which reduced the tensile stress in the slab just before it would have reached the tensile strength σ _t, thus resulting in a longer tensile length (Fig. 7b). The correlation length has a major influence on the strength of the weak layer, τ _p, with a shorter correlation length causing the peak to move closer to the bottom of the slope, resulting in a shorter tensile length than in the homogeneous case. However, cases with a longer correlation length (>25 m) have a similar behavior as cases with homogenous properties (Fig. 7a). The sharp transition in tensile length around 15–20 m of correlation length in our result (Fig. 8) is due to the sinus function itself, and this could be smoothed out by performing numerous simulations with sinus functions of different phases (i.e. by shifting the origin of function) or by using a GRF-based simulation (Appendix Fig. 12).

Figure 8. Sensitivity analysis for the tensile length L _t with regards to the mean $\overline {D}$, standard deviation S _D and correlation length $\epsilon$ of the slab depth sinus function.

To further investigate the results presented above in the case of a sinus function, we used GRF to simulate slab depth variation for a few GRF parameters. Figure 9 shows the distribution of tensile lengths from 50 realizations for four different sets of GRF parameters. The dashed line represents the homogeneous tensile length. For the small variance $S^2_D$ (0.005 m²), the distribution is narrow and centered around the homogeneous tensile length. For the higher $S^2_D$ (0.025 m²), the distribution of tensile length has a large variance, especially for the 5 m correlation length with values of tensile lengths ranging from 10 to 50 m. The distribution for a 1 m slab depth has nearly the same distribution shape, which is narrow for a small value of 0.005 m² $S^2_D$ and wider for a $S^2_D$ of 0.025 m² (Fig. 9). The variance of tensile length distribution is smaller for a 1 m mean slab depth compared to the 0.5 m distributions. The medians of tensile length distributions for the 1 m mean slab depth are larger compared to the homogeneous values denoted by the dashed line in (Fig. 9). The distributions are shifted towards bigger tensile length values compared to the homogeneous tensile length, with either a small or high variance (Fig. 9). This could be explained by a larger friction effect by the 1 m mean slab depth, causing a less important build-up of tensile stress in the slab which resulted in tensile length values larger compared to the homogeneous tensile length.

Figure 9. Probability density functions of tensile lengths L _t of 50 realizations for different mean slab depths $\overline {D}$ of 0.5–1 m, different variances $S^2_D$ of 0.005–0.025 m², and 5–15–30 m correlation lengths $\epsilon$. The dashed line represents the homogeneous tensile length for a mean slab depth of 0.5 m and 1 m.

Skier-triggering probabilities versus potential avalanche release sizes

This last result section presents both the sensitivity analysis of the skier-triggering probability and the potential avalanche release size. To obtain the potential size of the avalanche release, we multiplied the tensile lengths by the mean slab depth to get an estimate of the volume of snow mobilized for the avalanche. Typically the avalanche release size is computed from the volume of snow in movement but our simulations setup does not include the cross-slope length of the potential avalanche. Figure 10 shows both the probability of skier-triggering and the potential avalanche release size. The skier-triggering probability appears to be inversely related to the mean slab depth because the force induced by a skier at a given depth is inversely proportional to the depth. On the contrary, the potential avalanche release size increased with slab depth. We obviously got a correlation as we multiplied the tensile length by the mean slab depth to get the potential avalanche release size. However, the tensile strength in our simulations was parametrized based on the mean slab depth which also explains the longer tensile length values for thicker and stronger slabs. Furthermore, (Fig. 10) shows that the skier-triggering probability is increasing, with increasing the standard deviation for a given mean slab depth. However, the avalanche release size should be smaller when the standard deviation is increasing, but this result should be nuanced because (Fig. 9) shows that with a more natural spatial variability generated from a GRF, the avalanche size is quite variable compared to the homogeneous case but tends to be larger while the mean slab depth is increasing. The area in the heatmap where the maximum probability is correlated with the area of minimum potential avalanche release size, and the opposite is also present with the minimum of probability with the maximum size (Fig. 10). For the same mean slab depth values, high variation leads to a high skier-triggering probability but a lower potential avalanche release size. As the mean slab depth increases, the skier-triggering probability decreases to a point that for 0.9 and 1 m mean slab depth, the probability to trigger by a skier is almost impossible except with high slab depth variation, but the trigger could lead to relatively larger avalanches. These cases represent a scenario of low probability but high consequences with large avalanches, and this scenario is only possible with a slab depth spatial variation.

Figure 10. Probability of skier-triggering and normalized potential avalanche release size in relation to the mean slab depth and standard deviation of the slab depth. Potential avalanche release size combines the tensile length normalized with the largest tensile length multiplied by the mean slab depth. We show the standard deviation slab depth values for visual purposes but the variance values used with GRF method yield approximately the same values as the standard deviation slab depth used with the sinus function.

Discussion

Relevance of the study for practical implications

This study presents a novel mechanical–statistical approach to understanding the influence of slab depth spatial variability on the skier triggering probability of potential avalanche size. First, we showed a sensitivity analysis of the three parameters defining the variability. The increase of skier-triggering probability is inversely proportional to the mean slab depth, which was expected considering previous studies on this matter but without slab depth variation (Föhn, 1987; Monti and others, 2016; Gaume and Reuter, 2017). However, we show that the tensile length and the potential avalanche size are increasing proportionally to the mean slab depth, which is in agreement with previous studies showing that thicker and stronger slabs promote larger avalanches (Gaume and others, 2015). The variance and the correlation length results from the sensitivity analysis could be interpreted together. The combination of both creates some propagation spots (SPI below 1) on the fictional slope that skiers triggered. As the variance increases, it creates more propagation spots on the slope because it creates more areas where the slab is thinner so the skier can trigger the weak layer. A smaller correlation length leads to several small weak spots distributed across the slope compared to cases with longer correlation lengths, which led to fewer weak spots that cover more surface on the slope. The higher number of weak spots created using a smaller correlation length induced a higher skier-triggering probability than in the case with fewer weak spots covering more surfaces. In brief, a short-range variation creates numerous propagation spots leading to an increase in the skier triggering probability. Interestingly, several studies investigated the effect of the weak spots, created by heterogeneity in the weak layer, on natural release and obtained opposite conclusions (Gaume and others, 2013, 2014). The short-range variation of the weak layer creates several weak spots that are distributed across the slope like the short-range variation of slab depth. These short-range variations can be smoothed due to the elasticity of the slab (Gaume and others, 2013) and reduce the probability of a natural avalanche occurring. These studies also found that long-range variations create fewer weak spots but these weak spots were larger and covered more surfaces. This long-range variation was destabilizing the entire slope, in line with the so-called knock-down effect (Kronholm and Schweizer, 2003; Fyffe and Zaiser, 2004; Gaume and others, 2014). Finally, the variation of the snow mechanical properties (slab depth or weak layer cohesion), creates a spatial distribution of weak spots on the slope, and this variation does not have the same effect regarding the type of trigger between a skier (more unstable with a short-range variation) compared to a natural release (which is more unstable with a long-range variation).

This study also had an implicit objective to apply the method to test common knowledge that practitioners of the avalanche industry developed through many years of experience in the field. We realized that our method was heavily dependent on the skier style of the simulated skier trajectory. We show that a linear trajectory with a high down-slope radius R _down and a small cross-slope amplitude A _cross, reduce significantly the odds of triggering a weak spot. These results should be applied carefully in practice because it doesn't imply that skiing straight down the slope is ‘safer,’ whatever the term safer implicated in this context, but it only reduces the probability of triggering a weak spot resulting in an avalanche. The opposite trajectory was also simulated with a very high cross-slope amplitude and a small down-slope radius, representing an up-slope trajectory which we can translate to skinning up the up-slope with conversion (Fig. 5). Backcountry recreationists should not base their decision to ski a particular slope on their skiing trajectory. Decision-making in avalanche terrain is a complex task with many different aspects like terrain features, safety management, and other mountain hazards (Harvey and others, 2023)

We wanted to test if it was safer to ski closer to a preexisting ski track. Figure 6 shows a comparison between a structured approach that mimics skiing closely to a preexisting ski track from a completely random approach. The ECDF with a long correlation length had a median from 7 to 10 skiers on the slope before recording a trigger. It is important here to notice that in this method, the ECDF only represents cases where the first skier did not trigger, and then we started to add additional skiers both in a structured and random approach. The use of a safe first skier track is important because it mimics the fact that preexisting ski tracks could give information to other skiers that this trajectory did not trigger. Then, in spatial variation with a long correlation length like 30 m, the distance to the next weak spot could be on average 30 m away, explaining the fact that many skiers, with a spacing of 5 m could be on the slope before recording a trigger. These results confirm and quantify common knowledge used by ski guides and in avalanche awareness communication (Harvey and others, 2023).

Our results showed that the slab depth spatial variability adds randomness and unpredictability to skiing an avalanche-prone slope. A spatially homogeneous slab creates a binary outcome: either each skier triggers the avalanche or none do. But, the spatial variability of the slab depth creates a third regime in which some skiers trigger the avalanche and some do not on the same slope. The slab depth spatial variability creates weak spots on the slope which the skier trajectory will determine the outcome. The randomness is a result of the arbitrary anisotropic trajectory (down-slope) of the skier toward a potential weak spot on the slope. Our results also showed randomness in the tensile length obtained from our GRF simulations. The variation of slab depth sometimes induced an early crack arrest resulting in shorter tensile length, and sometimes the opposite outcome. The slab depth variation influences the crack speed during crack propagation. The observed increase in the speed of the crack caused a decrease in the tensile stress building up on the slab, which led to a longer tensile length. The same phenomenon was observed when the slab depth variation was slowing down the crack propagation, causing a sharp increase in tensile stress.

The last result presented in this study is the comparison between the skier triggering probability and the potential avalanche size. We show that it is more probable to trigger thinner and softer slabs compared to thicker and stronger slabs, but these thicker and stronger slabs could potentially create larger avalanche release sizes, also described and modeled by Gaume and others (2015). The increase of the slab depth standard deviation has an ambiguous effect because it increases the skier-triggering probability but reduces the avalanche release size. This latter was also observed by Gaume and others (2015) but regarding the variation of the weak layer cohesion instead of the slab depth on the avalanche release size. The relation between the likelihood to trigger by a skier compared to the propensity of propagation was initially described, through stability tests and field observations by van Herwijnen and Jamieson (2007). However, this study focused on fracture initiation and the propensity of propagation but not the potential avalanche size which refers to the dynamic crack propagation. Our statistical-mechanical provides a physically based validation and includes the effect of spatial variability of this well-known relationship (Fig. 10) which is the basis to describe the avalanche hazard in several countries (Statham and others, 2018; Techel and others, 2020).

Limits of the study and outlooks for future work

The novel methods used in this study present some limitations and assumptions that could be explored in future work. The analytical method considered the skier loading static along the skier tracks. In reality, the skier adds more pressure to the snow cover at the apex of each turn. This additional pressure in the apex could potentially trigger the weak layer in thicker slabs. Also, skier penetration was not taken into account, and this could potentially affect some of the presented results. The additional pressure coupled with the skier penetration depth could increase the skier-triggering probability using a skiing style with more turns compared to a linear skiing style. In addition, in the DAMPM model, the volumetric collapse of the weak layer and the induced slab bending is not taken into account because of the depth-averaged assumptions. In principle, at such a scale on an inclined slope, the slab tension represented by the DAMPM model is significantly higher than the slab bending which is anyway limited by the slab touchdown slab behind the crack tip (Benedetti and others, 2019).

We made the assumption that the weak layer strength would follow the slab depth variation locally. However, Bellaire and Schweizer (2011) has shown that is not always the case and the weak layer and the slab spatial pattern could differ. The relation between the slab depth and the weak layer strength is mainly due to the settlement and the friction from the slab weight, but some variation could still remain. The weak layer could be parameterized using the Mohr–Coulomb relation proposed by Gaume and others (2013), but with the friction term locally adjusted to the slab depth and the cohesion term could be set using a GRF with a different spatial pattern than the slab depth. This could lead to more realistic variations of snow properties and ultimately more realistic simulations. We choose not to follow this approach because we want to isolate only the slab depth variation to better understand individual drivers of the instability. Either approach should create weak spots triggered by skiers in the fictional slope, but maybe not in the same areas of the slope. However, the use of two different spatial patterns for the weak layer and the slab could create areas where the weak layer is stronger and the slab depth is thinner, the opposite of what this work is presenting. Such areas could promote slab tensile failure and potentially a crack arrest in the weak layer. Further investigations should be made on this matter to see if it could reduce the crack speed and potentially promote the crack to arrest, which we did not observe in our simulations. The covariance model used in the GRF did not contain any nugget. This choice created a smoother spatial realization of the slab depth without any ‘noise’ from the nugget. However, Kronholm and Birkeland (2005) has shown that increasing the nugget effect could promote crack arrest during dynamic crack propagation. This effect should be explored in future work regarding crack arrest and potential avalanche release size.

This study only investigates in numerical simulation the length at which the first tensile fracture occurred in the slab as a proxy to the potential avalanche release size. Tensile fracture is often related to crack arrest in the weak layer but not necessarily. There are no current studies to the best of our knowledge that explore the conditions of the crack arrest on an inclined snow slope. However, real-scale experiments on flat-terrain show some mechanical conditions when crack arrests were observed on 10 m long flat PST's (Bergfeld and others, 2021). They demonstrate some dissipation of energy during the dynamic crack propagation which could reduce the crack speed and possibly induce a crack arrest without a slab fracture. The dissipation of energy was due to the compression of the weak layer by the slab behind the crack tip. However, this phenomenon could be less important or may be absent on an inclined slope where the crack propagation changes to pure shear crack propagation (Trottet and others, 2022). Other dissipation mechanisms could influence the crack arrest like the softening behavior of the weak layer. Further investigation should focus on the effect of the softening on energy dissipation and crack arrest. Strong heterogeneity in the weak layer and topography changes should also be explored to further explain the conditions for a crack to arrest, and ultimately estimate the potential avalanche release size.

Conclusion

This study demonstrates the influence of slab depth spatial variability on the skier-triggering probability and the possible avalanche release area using a novel mechanical-statistical approach. We show that the spatial variability slab depth could increase the skier-triggering probability for thicker slabs when in a homogeneous case, a trigger by a skier is unlikely. For the possible avalanche release size, we show that the spatial variability of slab depth can induce a fluctuation in the tensile stress causing an early tensile failure resulting in smaller avalanches and also the opposite with bigger avalanches depending on the scale of the variability. We used the tensile length as a potential proxy for crack arrest but further research should focus on drivers for crack arrest dynamic crack propagation.

This study provides quantification and scientific evidence on the common knowledge that practitioners have developed throughout years of experience in the avalanche industry. We demonstrate the effect of skiing style on the probability to trigger an avalanche. We validate the idea with scientific evidence that skiing near a preexisting skier track could reduce your probability to trigger an avalanche compared to a random approach. This study demonstrates some processes during dynamic crack propagation regarding the variation of slab depth along a 1D slope. However, more research is needed to understand which drivers like topography or strong snow heterogeneity could potentially stop this dynamic crack propagation, both for the anticrack propagation in flat terrain and the super-shear regime on an inclined slope. Finally, this study shows, validates, and quantifies the well-known relationship between the likelihood and the size of an avalanche as well as common knowledge for safety guidelines in the avalanche community.