Snow fracture in relation to slab avalanche release: critical state for the onset of crack propagation

Abstract

The failure of a weak snow layer buried below cohesive slab layers is a necessary, but insufficient, condition for the release of a dry-snow slab avalanche. The size of the crack in the weak layer must also exceed a critical length to propagate across a slope. In contrast to pioneering shear-based approaches, recent developments account for weak layer collapse and allow for better explaining typical observations of remote triggering from low-angle terrain. However, these new models predict a critical length for crack propagation that is almost independent of slope angle, a rather surprising and counterintuitive result. Based on discrete element simulations we propose a new analytical expression for the critical crack length. This new model reconciles past approaches by considering for the first time the complex interplay between slab elasticity and the mechanical behavior of the weak layer including its structural collapse. The crack begins to propagate when the stress induced by slab loading and deformation at the crack tip exceeds the limit given by the failure envelope of the weak layer. The model can reproduce crack propagation on low-angle terrain and the decrease in critical length with increasing slope angle as modeled in numerical experiments. The good agreement of our new model with extensive field data and the ease of implementation in the snow cover model SNOWPACK opens a promising prospect for improving avalanche forecasting.

Introduction

Snow slab avalanches range among the most prominent natural hazards in snow-covered mountainous regions throughout the world. The winter 2014/15 served as a cruel reminder of the destructive power of this ubiquitous natural hazard with 132 fatalities, just for the European Alps. The ability to reliably forecast avalanche danger is therefore of vital importance and requires a sound understanding of avalanche release processes.

Avalanches are the result of numerous factors and processes interacting over a large range of temporal and spatial scales . While snow slab avalanches can come in many different sizes, from a few meters to several kilometers, they initiate within the snow cover by local damage processes at the grain scale. Indeed, the release of a dry-snow slab avalanche (Fig. a) requires the formation of a localized failure within a so-called weak layer (WL) buried below cohesive slab layers (Fig. b). The initial failure resulting in a macroscopic crack in the WL develops from micro-scale heterogeneities by damage accumulation or directly below a local overload such as a skier or a snowmobile . Stress concentrations at the crack tips will then determine whether crack propagation and eventually slope failure occurs , even when the average overlying stress is lower than the average weak layer strength knock-down effect;. The size of the initial crack at which rapid crack propagation occurs is called the critical crack length and represents an instability criterion for material failure . It is a crucial variable to evaluate snow slope instability .

(a) Crown fracture of a dry-snow slab avalanche in Mt. Baker, USA (^©Grant Gunderson). (b) Surface hoar weak layer. (c) Propagation saw test. The weak layer is represented in white, the slab in grey. The black dots are markers used for particle tracking to measure slab deformation.

[Figure omitted. See PDF]

Information on snow cover stratigraphy, especially the presence and characteristics of WLs and the overlying slab, is thus essential for avalanche forecasting. Traditionally, such information is obtained through manual snow cover observations, such as snow profiles and stability tests . However, these observations are time consuming, somewhat subjective, and only provide point observations. Snow cover models such as CROCUS and SNOWPACK provide a valuable alternative to obtain more highly resolved snow stratigraphy data. However, to evaluate snow slope instability based on model output, avalanche formation processes are greatly simplified and reduced to accounting for the balance between shear strength of the WL and shear stress due to the weight of the overlying slab, sometimes including a skier overload . This “strength-over-stress” approach is only relevant for failure initiation and does not account for crack propagation, the second fundamental process in avalanche release.

Due to the very complex nature of crack propagation in multilayered elastic systems under mixed-mode loading, theoretical and analytical approaches are not yet conceivable . In the past, simplifying assumptions have been used to propose analytical models for the critical crack length. For instance, , , and assumed a weak layer without thickness which allowed solution to the problem in the down-slope direction only, by neglecting the effect of the volumetric collapse of the weak layer as, for example, described by . In contrast, assumed a weak layer of finite thickness with a slope-independent failure criterion and a completely rigid behavior allowing to neglect the elastic mismatch between the slab and the weak layer. With the development of new field tests, in particular the propagation saw test (PST, Fig. c) , it is now possible to directly evaluate the critical crack length and thus determine crack propagation propensity. Particle tracking velocimetry (PTV) analysis of PSTs has highlighted the importance of the elastic bending of the slab induced by the loss of slab support due to weak layer failure (induced by a saw) prior to crack propagation . To include slab bending in the description of slab avalanche release mechanisms, proposed the anticrack model. This model provides an analytical framework to estimate the critical crack length as a function of slab properties (thickness, density, and elastic modulus) and the WL specific fracture energy, a WL property quantifying the resistance to crack propagation. While some crucial features of the mechanical behavior of the WL, including elasticity and shape of the failure envelope, are not included, the anticrack model provides a significant step forward as it accounts for various aspects that were left unexplained by previous theories, such as crack propagation on low-angle terrain and remote triggering of avalanches.

To evaluate the critical crack length based on the anticrack model, the WL specific fracture energy is required. Presently, it can be estimated using three existing methods: (i) through PTV or finite element analysis of the PST , (ii) from snow micro-penetrometer (SMP) measurements by integrating the penetration resistance over the thickness of the WL , and (iii) from X-ray computer tomography-based (CT) microstructural models . Depending on the method, estimates of the WL specific fracture energy can differ by as much as 2 orders of magnitude, resulting in widely different values of the critical crack length. Strength-of-material approaches have also been developed to evaluate the conditions for the onset of crack propagation . These methods require WL strength, a property which is more readily measurable , rather than the specific fracture energy. However, in contrast to the anticrack model, the latter strength-of-material approaches do not account for slab bending which leads to additional stress concentrations; hence these models tend to overestimate the critical length.

Clearly, the various methods to estimate the critical crack length all have their respective shortcomings, and a unified approach which incorporates all relevant processes is thus far missing. To overcome these limitations and take into account all the important physical ingredients, we propose to evaluate the critical crack length for different snowpack stratigraphies using discrete element simulations. Similar to the field experiments, in the simulations we gradually create a crack in the WL with a saw until rapid propagation occurs (Fig. ). On the basis of our numerical results, we then introduce a new expression for the critical crack length which accounts, for the first time, for the complex interplay between loading, elasticity, failure envelope of the WL, and its structural collapse. The predictive capabilities of this new expression, with respect to field data, are discussed and compared to previous models.

Successive snapshots (a–e) of a DEM simulation of the propagation saw test (PST). The plots on top of each snapshot represent illustrations of the shear stress $\mathit{\tau }$ (red line) in the WL. $D$ is the slab thickness (slope normal), $D_{\mathrm{wl}}$ is the WL thickness, $\mathit{\psi }$ is the slope angle, ${\mathit{\tau }}_{\max }$ is the maximum shear stress at the crack tip, ${\mathit{\tau }}_{\mathrm{p}}$ is the WL shear strength (dashed line), ${\mathit{\tau }}_{\mathrm{g}}=\mathit{\rho }gD\sin \mathit{\psi }$ is the shear stress due to the slab weight, and ${\mathit{\tau }}_{\mathrm{r}}$ is the residual frictional stress. $a_{\mathrm{c}}$ is the critical crack length, $\mathrm{\Lambda }$ is the characteristic length scale of the system, and $l_{\mathrm{0}}$ is the touchdown length (see Sect. 3). The red segment represents the saw used to cut inside the weak layer.

[Figure omitted. See PDF]

Methods

Discrete element model

We model crack propagation in a slab–WL system using the discrete element method (DEM). DEM is well suited to represent large deformations as well as the evolution of the microstructure of materials in a dynamic context . The simulations are performed using PFC2D (by Itasca), implementing the original soft-contact algorithm of . The numerical setup and the cohesive contact law implemented is fully described in . We recall here the main characteristics of the DEM model.

The simulated system (Fig. 2a) is 2-D and composed of a fixed substratum, a WL of thickness $D_{\mathrm{wl}}$ (varied between 0.02 and 0.06 $\mathrm{m}$), and a slab of thickness $D$ (varied between 0.2 and 0.8 $\mathrm{m}$). The slab is modeled with spherical elements of radius $r=0.01$ $\mathrm{m}$ with a square packing. As explained in , these elements are not intended to represent the real snow grains. They constitute entities of discretization used to model an elastic continuum of density $\mathit{\rho }$, Young's modulus $E$, and Poisson's ratio $\mathit{\nu }$. The WL is composed of elements of radius $r_{\mathrm{wl}}=r/\mathrm{2}$ with a packing of collapsible triangular shapes of the same size as the WL thickness (Fig. 2a) aimed at roughly representing the porous microstructure of persistent WLs such as surface hoar (Fig. 1b) or depth hoar.

We used the cohesive contact law detailed in . The bonds are characterized by specific elasticity and strength parameters which have been calibrated to obtain the desired macroscopic (bulk) properties. For the slab, numerical biaxial tests were performed to characterize the macroscopic Young's modulus $E$ as a function of micro-mechanical parameters. For the WL, mixed-mode shear-compression loading simulations were performed to determine the failure envelope (Fig. ). Through the triangular shape of the WL structure, the main features of real WL failure envelopes are captured, notably the closed envelope necessary to obtain failures both in shear and compression.

The applied loading represents a typical experimental setup of a PST . It consists of a combination of gravity (slope angle $\mathit{\psi }$) and advancing a rigid “saw” (in red in Fig. ) at a constant velocity $v_{\mathrm{saw}}=\mathrm{2}$ $\mathrm{m}{\mathrm{s}}^{-\mathrm{1}}$ through the WL. The saw thickness is $h_{\mathrm{saw}}=\mathrm{2}$ $\mathrm{mm}$ and the length of the system is $L=\mathrm{2}$ $\mathrm{m}$ .

Failure criterion FC${}_{\mathrm{1}}$ of our modeled weak layer (black circles) obtained from mixed-mode shear-compression loading tests. FC${}_{\mathrm{2}}$ is the high-rate mixed-mode failure envelope found by . The grey dotted lines represent angles of loading $\mathit{\psi }$ such as $\tan \mathit{\psi }={\mathit{\tau }}_{\mathrm{g}}/{\mathit{\sigma }}_{\mathrm{n}}$ where ${\mathit{\tau }}_{\mathrm{g}}$ is the shear stress. Compression corresponds to positive values of ${\mathit{\sigma }}_{\mathrm{n}}$.

[Figure omitted. See PDF]

Critical length $a_{\mathrm{c}}$ for crack propagation as a function of (a) Young's modulus $E$ of the slab, (b) slab density $\mathit{\rho }$, (c) slab thickness $D$, (d) WL thickness $D_{\mathrm{wl}}$, and (e) slope angle $\mathit{\psi }$. The symbols represent the critical length obtained from the DEM simulations and the solid lines represent the critical length modeled from Eq. () and for different failure behaviors. Dashed lines indicate the critical length obtained with the anticrack model assuming $w_{\mathrm{f}}=0.1$ $\mathrm{J}{\mathrm{m}}^{-\mathrm{2}}$.

[Figure omitted. See PDF]

Comparison with PST experiments

The data set consists of 93 PST experiments which were presented in . It includes the average slab density $\mathit{\rho }$, slab thickness $D$, slope angle $\mathit{\psi }$, and WL thickness $D_{\mathrm{wl}}$. The WL specific fracture energy $w_{\mathrm{f}}$ was evaluated from the penetration resistance of the weak layer using the SMP according to and ranges from 0.07 to 2.9 $\mathrm{J}{\mathrm{m}}^{-\mathrm{2}}$. showed a good correlation between the SMP-derived $w_{\mathrm{f}}$ and that derived using particle tracking and the anticrack model . The shear strength ${\mathit{\tau }}_{\mathrm{p}}$ of the WL was not measured but we used the mixed-mode shear-compression failure envelope defined by based on laboratory experiments. This failure envelope (in red in Fig. ), i.e., the relation between the shear strength ${\mathit{\tau }}_{\mathrm{p}}$ and the slope normal stress ${\mathit{\sigma }}_{\mathrm{n}}$, is described by the following Mohr–Coulomb–Cap model: $$\begin{array}{}& {\displaystyle }{\mathit{\tau }}_{\mathrm{p}}& {\displaystyle }={\mathit{\tau }}_{\mathrm{p}}^{\mathrm{mc}}=c+{\mathit{\sigma }}_{\mathrm{n}}\tan \mathit{\varphi }\mathrm{for}\mathit{\psi }>{\mathit{\psi }}_{\mathrm{t}},& {\displaystyle }{\mathit{\tau }}_{\mathrm{p}}& {\displaystyle }={\mathit{\tau }}_{\mathrm{p}}^{\mathrm{cap}}=b\sqrt{\mathrm{1}-{\displaystyle \frac{{\left({\mathit{\sigma }}_{\mathrm{n}},+,{\mathit{\sigma }}_{\mathrm{t}}\right)}^{\mathrm{2}}}{{\left({\mathit{\sigma }}_{\mathrm{c}},+,{\mathit{\sigma }}_{\mathrm{t}}\right)}^{\mathrm{2}}}}}\mathrm{for}\mathit{\psi }<{\mathit{\psi }}_{\mathrm{t}},\end{array}$$ where ${\mathit{\psi }}_{\mathrm{t}}=$ 23${}^{\circ }$ is the angle corresponding to a transition between the Mohr–Coulomb and the cap regimes, $c$ is the cohesion, $\mathit{\varphi }=$ 20${}^{\circ }$ is the friction angle, ${\mathit{\sigma }}_{\mathrm{t}}=c\tan \mathit{\varphi }$ is the tensile strength, ${\mathit{\sigma }}_{\mathrm{c}}=2.6$ $\mathrm{kPa}$ is the compressive strength, and $${\displaystyle }b=K\sqrt{{\displaystyle \frac{{\left({\mathit{\sigma }}_{\mathrm{t}},+,{\mathit{\sigma }}_{\mathrm{c}}\right)}^{\mathrm{2}}}{{\left({\mathit{\sigma }}_{\mathrm{t}},+,{\mathit{\sigma }}_{\mathrm{c}}\right)}^{\mathrm{2}}-{\left(\frac{K}{\tan \mathit{\varphi }}\right)}^{\mathrm{2}}}}}.$$ $K=\mathrm{1}$ $\mathrm{kPa}$ is the maximum shear strength . The cohesion $c$ (shear strength for ${\mathit{\sigma }}_{\mathrm{n}}=\mathrm{0}$) can be derived from the WL specific fracture energy $w_{\mathrm{f}}$ using the results of : $${\displaystyle }c={\displaystyle \frac{\sqrt{\mathrm{2}D{E}^{\prime }{w}_{\mathrm{f}}}}{\mathrm{2}\mathrm{\Lambda }}},$$ where $\mathrm{\Lambda }$ is a characteristic length scale of the system see Sect. 3 and. Note that, for the 93 PST experiments, the slope normal stress ${\mathit{\sigma }}_{\mathrm{n}}$ was lower than 2 $\mathrm{kPa}$ and thus only the Mohr–Coulomb part of the failure envelope (Eq. 1) was used to compute the shear strength ${\mathit{\tau }}_{\mathrm{p}}$. For these stress states (low slope normal stress), showed that failure was almost independent of the loading rate (in the brittle range) and thus independent of fast sintering effects .

The Young's modulus of the slab $E$, which was not measured, was derived from density according to : $${\displaystyle }E=5.07\times {10}^{\mathrm{9}}{\left({\displaystyle \frac{\mathit{\rho }}{{\mathit{\rho }}_{\mathrm{ice}}}}\right)}^{5.13},$$ with ${\mathit{\rho }}_{\mathrm{ice}}=917$ $\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$. The WL shear modulus $G_{\mathrm{wl}}$ was taken constant equal to 0.2 $\mathrm{MPa}$ according to the laboratory experiments performed on snow failure by and Poisson's ratio of the slab $\mathit{\nu }$ was taken equal to 0.2 .

Results

DEM simulations

In the simulations, the crack of length $a$ created by the advancing saw in the WL induces slope-parallel and slope-normal displacements of the slab due to gravity leading to tension and bending in the slab. This results in stress concentrations at the crack tip where the shear stress $\mathit{\tau }={\mathit{\tau }}_{\max }$ is maximum and larger than the shear stress due to slab weight ${\mathit{\tau }}_{\mathrm{g}}$. The critical crack length $a_{\mathrm{c}}$ required for the onset of dynamic crack propagation in the WL is reached when ${\mathit{\tau }}_{\max }$ meets the shear strength ${\mathit{\tau }}_{\mathrm{p}}$ (Fig. c).

We performed a series of systematic simulations to investigate the influence of snow cover parameters on $a_{\mathrm{c}}$ (Fig. ). Slab properties (slab density $\mathit{\rho }$, slab elastic modulus $E$, slab thickness $D$), WL thickness $D_{\mathrm{wl}}$, and slope angle $\mathit{\psi }$ were varied independently in the simulations. Overall, $a_{\mathrm{c}}$ was found to increase with increasing elastic modulus of the slab $E$ and with WL thickness $D_{\mathrm{wl}}$. On the contrary, $a_{\mathrm{c}}$ decreased with increasing slab density $\mathit{\rho }$, with increasing slab thickness $D$ and with increasing slope angle $\mathit{\psi }$.

Analytical expression for the critical crack length

The discrete element simulations revealed that the maximum shear stress at the crack tip can be decomposed into two terms related to slab tension (${\mathit{\tau }}_{\max }^{\mathrm{t}}$) and slab bending (${\mathit{\tau }}_{\max }^{\mathrm{b}}$): $${\displaystyle }{\mathit{\tau }}_{\max }={\mathit{\tau }}_{\max }^{\mathrm{t}}+{\mathit{\tau }}_{\max }^{\mathrm{b}}.$$ When disregarding slab bending (weak layer with no thickness), the maximum stress ${\mathit{\tau }}_{\max }^{\mathrm{t}}$ depends on the shear stress due to the weight of the slab ${\mathit{\tau }}_{\mathrm{g}}$, the crack length $a$, and a characteristic length scale of the system $\mathrm{\Lambda }$ : $${\displaystyle }{\mathit{\tau }}_{\max }^{\mathrm{t}}={\mathit{\tau }}_{\mathrm{g}}\left(\mathrm{1},+,{\displaystyle \frac{a}{\mathrm{\Lambda }}}\right).$$ The length scale $\mathrm{\Lambda }$ represents the characteristic scale of the exponential decay of the shear stress $\mathit{\tau }$ close to the crack tip (Fig. b). It is given by $\mathrm{\Lambda }={\left(E^{\prime },D,D_{\mathrm{wl}},/,G_{\mathrm{wl}}\right)}^{\mathrm{1}/\mathrm{2}}$, where $E^{\prime }=E/(\mathrm{1}-{\mathit{\nu }}^{\mathrm{2}})$ is the plane stress elastic modulus of the slab and $G_{\mathrm{wl}}$ the WL shear modulus . We assume the shear stress inside the WL to be equal to the gravitational stress acting at the slab–WL interface, i.e., ${\mathit{\tau }}_{\mathrm{g}}=\mathit{\rho }gD\sin \mathit{\psi }$. Note that in the limiting case of a WL with zero thickness ($D_{\mathrm{wl}}\to \mathrm{0}$), the characteristic length is defined as $\mathrm{\Lambda }={\left(E^{\prime },D,/,k_{\mathrm{wl}}\right)}^{\mathrm{1}/\mathrm{2}}$, with $k_{\mathrm{wl}}$ the shear stiffness of the interfacial WL. Hence, as in the anticrack model (where WL failure is considered as an interfacial failure), WL thickness $D_{\mathrm{wl}}$ plays no role in the model for a constant WL stiffness $k_{\mathrm{wl}}$.

The tension term alone is unable to predict stress concentrations and thus crack propagation on flat terrain ($\mathit{\psi }=\mathrm{0}$), a process that exists, exemplified by numerous field observations and our DEM simulations (Fig. e). To resolve this discrepancy, the second term in Eq. () accounts for slab bending induced by WL collapse. Our DEM simulations showed that this term depends on the slope normal stress ${\mathit{\sigma }}_{\mathrm{n}}$ and the ratio $a/\mathrm{\Lambda }$ (Fig. b) and can be expressed as $${\displaystyle }{\mathit{\tau }}_{\max }^{\mathrm{b}}\approx {\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\mathit{\sigma }}_{\mathrm{n}}{\left({\displaystyle \frac{a}{\mathrm{\Lambda }}}\right)}^{\mathrm{2}}.$$ For the comparison with the anticrack model which assumes a rigid weak layer, one can consider the bending of a beam over a rigid foundation . In this case, the bending term ${\mathit{\tau }}_{\max }^{\mathrm{b}}$ would scale with ${\mathit{\sigma }}_{\mathrm{n}}(a/D{)}^{\mathrm{2}}$, independent of the elastic properties of the slab and the WL (similar to the anticrack model). In the present formulation, scaling with $a/\mathrm{\Lambda }$ instead of $a/D$ provides a means to account for the elastic mismatch between the slab and the WL and to adequately reproduce the numerical results (Fig. ).

Ratio between the shear strength ${\mathit{\tau }}_{\mathrm{p}}$ and the slope normal stress ${\mathit{\sigma }}_{\mathrm{n}}$ vs. the ratio between the critical length $a_{\mathrm{c}}$ and (a) slab thickness $D$ or (b) characteristic length $\mathrm{\Lambda }$ for flat terrain ($\mathit{\psi }={\mathrm{0}}^{\circ }$, i.e., ${\mathit{\tau }}_{\mathrm{g}}=\mathrm{0}$). The symbol/color in the legend indicates the parameter which was varied in the DEM simulations. The dashed line corresponds to Eq. ().

[Figure omitted. See PDF]

From Eq. () the critical length can be obtained by solving ${\mathit{\tau }}_{\max }={\mathit{\tau }}_{\mathrm{p}}$, where ${\mathit{\tau }}_{\mathrm{p}}$ is the shear strength given by the failure envelope of the material : $${\displaystyle }{a}_{\mathrm{c}}=\mathrm{\Lambda }\left[{\displaystyle \frac{-{\mathit{\tau }}_{\mathrm{g}}+\sqrt{{\mathit{\tau }}_{\mathrm{g}}^{\mathrm{2}}+\mathrm{2}{\mathit{\sigma }}_{\mathrm{n}}\left({\mathit{\tau }}_{\mathrm{p}},-,{\mathit{\tau }}_{\mathrm{g}}\right)}}{{\mathit{\sigma }}_{\mathrm{n}}}}\right].$$ Theoretically, this expression is valid only when crack propagation occurs before the slab touches the broken WL, i.e., when the vertical displacement induced by bending remains lower than the collapse height $h_{\mathrm{c}}$. The length $l_{\mathrm{0}}$ (Fig. d) required for the slab to come into contact with the broken WL can be expressed using beam theory: $l_{\mathrm{0}}={\left(\frac{\mathrm{2}E{D}^{\mathrm{2}}{h}_{\mathrm{c}}}{\mathrm{3}\mathit{\rho }g\cos \mathit{\psi }}\right)}^{\mathrm{1}/\mathrm{4}}$ . For realistic model parameters, $a_{\mathrm{c}}$ was always substantially lower than $l_{\mathrm{0}}$ (not shown).

The agreement between Eq. () and results from the DEM simulations is excellent (red solid lines in Fig. ). We emphasize that scaling of ${\mathit{\tau }}_{\max }^{\mathrm{b}}$ with $a/\mathrm{\Lambda }$ is of critical importance. It also provides an explanation for the gentler decrease of $a_{\mathrm{c}}$ with $D$ compared to $\mathit{\rho }$, even though $D$ and $\mathit{\rho }$ equally contribute to the load. Indeed, for a constant load, thicker slabs will result in lower stress concentrations at the crack tip (Eq. ) due to an increase of $\mathrm{\Lambda }$.

Comparison between measured and modeled critical crack lengths using the anticrack model (black circles) and our new model (Eq. , red stars). The continuous lines represent linear fits, both highly significant with $p<$ 0.01; in black, $a_{\mathrm{c}}={\mathit{\gamma }}_{\mathrm{H}}{a^{\prime }}_{\mathrm{c}}+{\mathit{\delta }}_{\mathrm{H}}$ with ${\mathit{\gamma }}_{\mathrm{H}}=0.46$ (0.32, 0.61), ${\mathit{\delta }}_{\mathrm{H}}=0.23$ (0.21, 0.30), and $R_{\mathrm{H}}^{\mathrm{2}}=0.24$; in red, $a_{\mathrm{c}}={\mathit{\gamma }}_{\mathrm{G}}{a^{\prime }}_{\mathrm{c}}+{\mathit{\delta }}_{\mathrm{G}}$ with ${\mathit{\gamma }}_{\mathrm{G}}=1.023$ (0.872, 1.173), ${\mathit{\delta }}_{\mathrm{G}}=0.03$ ($-$0.01, 0.07), and $R_{\mathrm{G}}^{\mathrm{2}}=0.58$. The numbers in brackets represent the 95 % confidence interval. The dashed line represents the $\mathrm{1}:\mathrm{1}$ line.

[Figure omitted. See PDF]

The predictions of Eq. () also compare well with results obtained from 93 PST experiments (Fig. ). Overall, our model provides very good estimates of the measured critical crack lengths, as demonstrated by the proximity of the data to the $\mathrm{1}:\mathrm{1}$ line despite substantial scatter ($R^{\mathrm{2}}=0.58$). As for the simulations, the critical length in PSTs was always lower than the length $l_{\mathrm{0}}$ (not shown).

Discussion

Comparison with the anticrack model

We compare how well our new analytical expression (Eq. ) and the anticrack model can reproduce the dependence of the critical crack length on system properties as obtained with our DEM simulations (Fig. ). The anticrack model reproduces the influence of $E$, $\mathit{\rho }$, and $D$ on $a_{\mathrm{c}}$ well for $\mathit{\psi }=\mathrm{0}$, although less accurately than Eq. (). However, the influence of WL thickness $D_{\mathrm{wl}}$ and slope angle $\mathit{\psi }$ on $a_{\mathrm{c}}$ was very poorly reproduced by the anticrack model, both in terms of absolute values and trends. In particular, a slope angle $\mathit{\psi }>\mathrm{0}$ would lead to similar trends of $a_{\mathrm{c}}$ with $E$, $\mathit{\rho }$ and $D$ but with overestimated values.

The decrease of $a_{\mathrm{c}}$ with slope angle, observed in our DEM results and predicted by Eq. (), is of particular interest. This trend is in clear contradiction with one of the main outcomes of the anticrack model , namely that the critical length is almost independent of slope angle. The discrepancy arises from the fact that the anticrack model (i) assumes that the failure behavior of the WL is slope independent, (ii) disregards WL elasticity, and (iii) does not adequately account for the interplay between tension and bending in the slab as also shown in . Concerning WL thickness, a thin WL leads to higher stress concentrations in bonds between the grains and thus to a smaller critical crack length (Fig. d). This effect cannot be reproduced by the anticrack model due to the rigid character of the WL.

For low-angle terrain, the anticrack model and our new formulation yield similar results. However, this is where the similarities end. Indeed, overall the anticrack model overestimates $a_{\mathrm{c}}$ and more closely resembles a model which only accounts for stresses due to slab bending: $a_{\mathrm{c}}^{\mathrm{b}}=\mathrm{\Lambda }\sqrt{{\mathit{\tau }}_{\mathrm{p}}/{\mathit{\sigma }}_{\mathrm{n}}}$ (obtained by solving ${\mathit{\tau }}_{\max }^{\mathrm{b}}={\mathit{\tau }}_{\mathrm{p}}$). For steep slopes ($\mathit{\psi }>$ 30${}^{\circ }$), where the shear stress at the crack tip due to slab bending becomes negligible compared to that due to slab tension, critical crack length values obtained from Eq. () strongly differ from the prediction of the anticrack model and converge on the contrary towards a purely tensile model, generally referred to as “pure shear model”: $a_{\mathrm{c}}^{\mathrm{t}}=\mathrm{\Lambda }({\mathit{\tau }}_{\mathrm{p}}/{\mathit{\tau }}_{\mathrm{g}}-\mathrm{1})$ (obtained by solving ${\mathit{\tau }}_{\max }^{\mathrm{t}}={\mathit{\tau }}_{\mathrm{p}}$, Fig. e).

illustrated the low dependence of the critical crack length on slope angle with results from field experiments presented by . However, these PST experiments were performed on a non-persistent WL consisting of precipitation particles and measurements made on the flat were performed 1 day before the experiments made on slopes . This indicates that the trend with slope angle may be influenced by the burial time of the WL since sintering and settlement effects can strongly affect snowpack properties within 1 day, especially with the layer of precipitation particles which was tested . Furthermore, assumed snow cover properties independent of slope angle, which is somewhat questionable since snowpack properties can also change with slope angle, thus obscuring the true slope angle influence. As an example, for their validation, assumed a constant slab thickness $D=11$ $\mathrm{cm}$ over the different slope angles $\mathit{\psi }$, while $D$ decreased with increasing $\mathit{\psi }$ according to . In addition, it is also known that weak layer strength , slab density , and thus the elastic modulus are strongly depend on slope angle. Hence we argue that the dependence of the critical crack length on slope angle obtained from a model with fixed value of the other parameters should not be compared to the trend observed in the experiments, which is the result of a combination of many varying properties. Instead, one should directly compare the measured critical crack length to the modeled one, taking as input parameters the properties measured at the location where the PST was performed.

By comparing the anticrack model to the 93 PST measurements (Fig. ), we see that $a_{\mathrm{c}}$ is generally overestimated, especially for short critical crack lengths and steep slopes (${35}^{\circ }<\mathit{\psi }<{45}^{\circ }$). For higher values of $a_{\mathrm{c}}$ and gentler slopes, the anticrack predictions better agree with our formulation, even though they still remain mostly above the $\mathrm{1}:\mathrm{1}$ line.

Slope angle dependence

We showed that the critical crack length $a_{\mathrm{c}}$ decreases with increasing slope angle $\mathit{\psi }$ for a PST with slope-normal faces, a constant slab thickness $D$, and constant values of the mechanical properties. However, the rate of decrease of $a_{\mathrm{c}}$ with $\mathit{\psi }$ is strongly influenced by the elastic modulus $E$ and thickness $D$ of the slab. Low values of $E$ and/or $D$ lead to a gentler decrease of $a_{\mathrm{c}}$ with $\mathit{\psi }$ (Fig. ).

Effect of the slab elastic modulus on the slope angle dependency of the critical crack length (Eq. ) for $\mathit{\rho }=200$ $\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$, $D=0.2$ $\mathrm{m}$, $D_{\mathrm{wl}}=\mathrm{4}$ $\mathrm{cm}$. Inset: effect of slab thickness on the slope angle dependency of the critical crack length for $\mathit{\rho }=200$ $\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$, $E=\mathrm{2}$ $\mathrm{MPa}$, $D_{\mathrm{wl}}=\mathrm{4}$ $\mathrm{cm}$.

[Figure omitted. See PDF]

(a) Critical crack length vs. slope angle: comparison between the data of (black circles) and our new model (Eq. , red line), the anticrack model (purple dashed-line), and a pure shear model (Eq. , green dotted line) for a constant slab depth $H=0.35$ $\mathrm{m}$ ($D=H\cos \mathit{\psi }$) and the same input parameters as in with a semi-log scale. The cohesion $c=500$ $\mathrm{Pa}$ was estimated based on the hand hardness index provided in using the parametrization by and . Inset: linear scale. (b) Effect of geometry on the slope angle dependency for the PST. SNF: slope normal faces; VF: vertical faces; const.: constant.

[Figure omitted. See PDF]

However, if slab depth $H$ (vertical) is constant with respect to slope angle, the slab thickness decreases with increasing slope angle according to $D=H\cos \mathit{\psi }$. Since a lower slab thickness leads to a higher critical crack length (Fig. 4c), this effect leads to an apparent reduction of the decrease of $a_{\mathrm{c}}$ with $\mathit{\psi }$. As an illustration, we compare our model (Eq. ) to the PST experiments presented for which the slab density and elastic modulus were very low (storm snow, $\mathit{\rho }=84$ $\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}$, $E=0.22$ $\mathrm{MPa}$; Fig. ). The low elastic modulus thus leads (Eq. ) to a very gentle decrease of $a_{\mathrm{c}}$ with $\mathit{\psi }$ in line with the experimental data. The anticrack model was also plotted in Fig. a and shows very comparable results. However, the values of the WL specific fracture energy $w_{\mathrm{f}}$ and slab elastic modulus $E$ in were estimated by a fit of the anticrack model to the data using the method described by which explains the good agreement. Interestingly, also a pure shear model (Eq. ) with the same input parameters as for our model (Eq. ) would lead to a reasonable agreement for steep slopes ($\mathit{\psi }>$ 30${}^{\circ }$). In the studies of and , the significant difference obtained between the anticrack model and the pure shear model is an artifact simply due to the fact that the same specific fracture energy was taken as input for both models although the underlying physical assumptions are strictly incompatible: the pure shear model considers a quasi-brittle behavior for the weak layer and the anticrack model considers a purely rigid behavior. In fact, for $\mathit{\psi }>$ 30${}^{\circ }$ and short critical crack lengths which are typically encountered in field experiments, recently showed from the energy balance equations that both approaches lead to very comparable results, which is confirmed by our new model.

Finally, geometrical effects significantly influence how the critical crack length depends on slope angle. Figure b shows the critical crack length as a function of slope angle for three different PST configurations: (i) constant slab thickness $D$ and slope normal faces (SNF); (ii) constant slab depth $H$ and SNF; (iii) constant slab depth and vertical faces (VF). The vertical character can be accounted for by adding $D/\mathrm{2}\tan \mathit{\psi }$ to the critical crack length as proposed by (see Supplement of ). We clearly observe that the decrease of $a_{\mathrm{c}}$ with $\mathit{\psi }$ is gentler with a constant slab depth $H$ than with a constant slab thickness $D$ as shown before. In addition, we observe an increase of the critical crack length with increasing slope angle when the PST is made with VF and when the slab depth is constant. This is in line with the PST experiments of performed with VF and a constant slab depth $H$. It seems that accounted for neither the vertical character of the faces nor the decrease of slab thickness with slope angle in their comparison to the data of . The increasing trend predicted by our model with a constant slab depth $H$ and VF might explain why the extended column test (ECT) scores were observed to increase with increasing slope angle .

Relevance and limitations

Performing DEM simulations allowed us to investigate crack propagation in weak snow layers without relying on the same strong assumptions concerning the weak layer as previous research . For the sake of developing theoretical models, these studies considered either a purely interfacial weak layer or a weak layer composed of a completely rigid material with a slope-independent failure criterion . On the contrary, in our simulations, the weak layer is characterized by a finite thickness, an elasticity, and a mixed-mode failure envelope in line with results of recent laboratory experiments . These DEM simulations can thus be seen as numerical laboratory experiments in which the effect of slab and weak layer properties on crack propagation can be investigated independently (which is impossible to do in the field) and from which analytical expressions can be inferred using a strength-of-material approach. This important step forward allows us to reconcile the shear- and collapse-based approaches. For example, our model can describe crack propagation in flat terrain providing the same results as the anticrack model. Furthermore, it predicts the decrease of the critical crack length with increasing slope angle in line with shear-based models and in contrast with the anticrack model since the latter assumes rigidity and slope-independent failure of the weak layer. Note that in the simulations and in reality, slab bending also induces shear stresses within the slab, leading to possible slope normal stress variations in the WL. This effect is not accounted for in our analysis. However, the good agreement between Eq. () and DEM results (Fig. ) suggests that it is in fact of second order, thereby validating the assumption that the maximum shear stress at the crack tip has two main contributions related to slab tension and bending (Eq. ).

(a) Seasonal profile of the simulated critical crack length (winter 2014–2015) at Steintälli (Davos, Switzerland) on the flat. (b) Vertical profile of the critical crack length (modeled and from field PSTs) and SSI/10 for the date marked by the vertical red line in (a). The grain type is shown on the right following .

[Figure omitted. See PDF]

In a recent study, showed that the DEM model can also reproduce the dynamic phase of crack propagation as well as fracture arrest in the slab which was treated as an elastic–brittle material. In particular, the crack propagation speed and distances obtained by PTV analysis of the PST were well reproduced. It was also shown that the propagation distance (distance between the lower edge and slab fracture) was almost always higher than the critical crack length except for combinations of very low slab densities and thicknesses. This behavior is also observed in field experiments. Accordingly, treating the slab as a linear elastic material before the onset of crack propagation is justified. This assumption was also confirmed by recent field studies showing that the slab displacement obtained with particle tracking can be described by beam theory with a linear elastic assumption. Hence, with the present study, we show that our DEM model is able to address the whole crack propagation process.

The main limitation of our model is the uniform character of the slab. In this paper, the multilayered character of the slab was not accounted for, for clarity reasons since the phenomenon is already very complex. However, the elastic moduli of the slab layers have a very important influence on slab deformation and thus on the critical crack length . For the comparison with the experiments, the elastic modulus was computed from the average slab density. However, in practice, a slab with a uniform density $\mathit{\rho }$ will deform differently than a slab of average density $\mathit{\rho }$ consisting of various layers with contrasting properties. This is probably the reason why significant scattering is observed in Fig.  although the overall agreement is good.

Concerning the weak layer, the schematic microstructure considered in this study is sufficient to capture the main features of the failure envelope . Considering more complex microstructures for the weak layer might lead to a better quantitative agreement with experimental data. For instance, it was shown that with a weak layer produced by ballistic deposition, the experimental failure envelope obtained by could be reproduced. In the future, performing numerical simulations accounting for the real microstructure of weak snow layers, as derived from X-ray microtomographic images , represents an interesting prospect. Nevertheless, if such refinements can certainly have an impact on the shear strength value ${\mathit{\tau }}_{\mathrm{p}}$, they should not change the structure of Eq. ().

Another important aspect is the relevance of our new model with regards to slab avalanche release. We showed that our model was able to reproduce crack propagation at the scale of the PST. However, at the slope scale, 3-D effects, slope-transverse propagation, terrain, and snowpack variability might make the process even more complex. Nevertheless, it was shown that the critical crack length correlates very well with signs of instability . In particular, they showed that no signs of instability were recorded for $a_{\mathrm{c}}>0.4$ $\mathrm{m}$ while whumpfs, cracks, and avalanches were observed for $a_{\mathrm{c}}<0.4$ $\mathrm{m}$. Hence, our new model of critical crack length can be of major importance in view of avalanche forecasting.

Application to simulated snow stratigraphy

The snow cover model SNOWPACK , which simulates the temporal evolution of snow stratigraphy, is used for operational avalanche forecasting in Switzerland. Potential weak layers in the simulated snow profiles are identified by calculating the structural stability index (SSI), an index based on the balance between shear stress and shear strength . The SNOWPACK model also provides all necessary variables to determine the critical crack length based on Eq. (). To demonstrate the practical applicability, we performed a simulation for the 2014–2015 winter at the location of an automatic weather station above Davos, Switzerland (Fig. ). Note that the critical length was arbitrarily set to 1 $\mathrm{m}$ in the first 10 $\mathrm{cm}$, since avalanche probability for such shallow layers is generally very low . The same was done when computed values of the critical length exceeded 1 $\mathrm{m}$. Short critical crack lengths clearly highlight potential WLs in the snowpack during the season (Fig. a). At the end of the dry-snow season, around 10 April, the percolation of liquid water into the snow cover resulted in a rapid increase in shear strength and thus in larger critical crack lengths throughout the snow cover.

On 3 March 2015 we performed several PSTs on three WLs at the location of the automatic weather station. The SNOWPACK simulation for that specific day clearly shows local minima in the calculated critical crack length for these three WLs (Fig. b). Modeled critical crack lengths were in good agreement with PST field measurements (black circles in Fig. b), and SNOWPACK was able to reproduce the observed increase in $a_{\mathrm{c}}$ with increasing depth of the WL. recently followed the temporal evolution of the critical cut length and showed that the implementation of Eq. () is very sensitive to the parametrization of ${\mathit{\tau }}_{\mathrm{p}}$ used in SNOWPACK . Finally, layers for which critical crack lengths were lower generally also corresponded to layers with local minima in the SSI, suggesting that a combination of SSI and $a_{\mathrm{c}}$ may provide a more reliable instability criterion .

Conclusions

We proposed a new analytical expression to assess the conditions for the onset of crack propagation in weak snowpack layers. The formulation was developed based on discrete element simulations; it accounts for crucial physical processes involved in crack propagation in snow, namely the complex mechanical behavior of the WL and the mixed stress states in the slab induced by slab tension and bending resulting from WL collapse. A critical parameter in the formulation is the length scale $\mathrm{\Lambda }$, which accounts for the elastic mismatch between the slab and the WL.

The analytical expression for the critical crack length reproduced field data obtained with 93 PST experiments. In contrast, the anticrack model , although appropriate for flat terrain, significantly overestimated the critical length for steep slopes, where avalanches release. Furthermore, our model predicts that the critical crack length decreases with increasing slope angle. This shows that triggering an initial failure leading to slab avalanche release is more likely on steep rather than on low-angle slopes, a rather intuitive result. Nevertheless, our model still allows for crack propagation on flat terrain and remote triggering of avalanches, both of which are widely documented by countless field observations.

Finally, our new expression was implemented in the snow cover model SNOWPACK to evaluate the critical crack length for all snow layers throughout the entire season. While validation is still required, this opens promising perspectives to improve avalanche forecasting by combining traditional stability indices with a new metric to evaluate crack propagation propensity.

Data availability

The critical crack length model is implemented in the SNOWPACK model, which is available under the GNU Lesser General Public Licence Version 3 and can be retrieved at http://models.slf.ch. PST data are available upon request from Johan Gaume ([email protected]).

The authors declare that they have no conflict of interest.

Acknowledgements

We are grateful to all SLF colleagues who assisted in field data collection. We thank Benjamin Reuter for the SMP-derived specific fracture energy data and for insightful discussions and comments on the paper. We acknowledge the constructive comments of two anonymous reviewers as well as Ned Bair who helped us to improve our paper. Johan Gaume has been supported by the Ambizione grant of the Swiss National Science Foundation (PZ00P2_161329). Edited by: E. Larour Reviewed by: E. H. Bair and two anonymous referees