Large-strain elastoplastic model

To model the observed process of anticrack propagation in snow, we developed a large-strain elastoplastic model. Material deformation is characterized by the strain measure. Assuming there is a deformation map ϕ(X, t) that maps undeformed coordinate X to a deformed coordinate x, the deformation gradient F is defined as ∂ϕ/∂X. Our physical model assumes finite strain elastoplasticity, where F is decomposed into elastic (FE) and plastic (FP) parts as F = FEFP (multiplicative elastoplasticity). The elastic deformation gradient is computed using the isotropic Hooke’s law of elasticity (see Methods section for more details).

For plasticity, the yield function y(τ) ≤ 0 defines admissible stress states in an elastoplastic continuum. We model snow based on the critical state plasticity theory for soil mechanics25,26. For any stress τ, there exist a mean effective stress (or pressure) p and a deviatoric stress s. They are given by

$$p = - \frac{1}{d}{\mathrm{tr}}({\boldsymbol{\tau }}),$$ (1)

$${\bf{s}} = {\boldsymbol{\tau }} + p{\bf{I}},$$ (2)

respectively, where d = 2 or 3 is the problem dimension, I is the identity matrix and compression corresponds to p > 0. According to the Von Mises theory27, we can derive the Mises equivalent stress q, given by q = (3/2 s : s)1/2 (so that \(q = \left| {\tau _1 - \tau _2} \right|\) for 2D and \(q\) = \(\sqrt {\frac{1}{2}\left( {(\tau _1 - \tau _2)^2 + (\tau _2 - \tau _3)^2 + (\tau _3 - \tau _1)^2} \right)}\) for 3D, in principal stress space).

Recent experiments11 and simulations based on X-ray microtomography28,29,30,31 highlighted the mixed-mode nature of snow failure including tensile, shear and compression failure modes. Given these past studies, it appears that an ellipsoid yield function is appropriate to reproduce this mixed-mode character. Hence we chose to start from the modified cam clay (MCC) yield surface32 which has been widely used in the area of soil mechanics. Note that the analogy between snow and clay was already made by McClung13 who extended the clay model of Palmer and Rice33 to model shear fractures induced by strain softening. However, the MCC model is originally cohesionless and does not exhibit any stress under extension, similar to dry sand. Hence, cohesion was added to the yield function by shifting the MCC model along the p-axis. We thus propose a new cohesive cam clay (CCC) model similar to that of Meschke et al.34 with the following yield surface:

$$y(p,q) = q^2(1 + 2\beta ) + M^2\left( {p + \beta p_0} \right)\left( {p - p_0} \right),$$ (3)

where p 0 represents the consolidation pressure and M is the slope of the cohesionless critical state line (CSL) that controls the amount of friction, β represents the ratio between tensile and compressive strength and controls the amount of cohesion (β ≥ 0). This yield surface is represented in Fig. 1a. Both MCC and our model are ellipsoids and are symmetric around the hydrostatic axis.

Fig. 1 Overview of the elastoplastic model. a Cohesive (black line) and cohesionless (dashed gray line) cam clay yield surface in the p–q space. The red line corresponds to the Critical State Line. b Illustration of the hardening models \(p_0\left( {\epsilon _V^P} \right)\) (for the slab) and p 0 (η) (for the weak layer): the black arrow shows the classical hardening law used for the snow slab in which p 0 increases in compression \(\left( {\dot \epsilon _V^P < 0} \right)\); the blue arrows represent the new softening model for the weak layer for which p 0 decreases under compression \(\left( {\dot \eta = \alpha \left| {\dot \epsilon _V^P} \right| > 0} \right)\) until \(\epsilon _V^P = 0\) after which the classical hardening law is used with β = 0. c Typical p–\(\epsilon _V\) curve obtained for the unconfined compression of the weak layer in experiment number 2 (see Methods section for model parameters) for the classical hardening law (in black) and the new softening one (in blue). d Same as c but for the q–\(\epsilon _V\) curve. In c and d, p and q in the weak layer (blue curves) do not perfectly reach zero after softening due to a loss of homogeneity (failure localization) Full size image

For the dense snow slab, the hardening and softening is modeled by expanding and shrinking the yield surface which is performed by varying p 0 . We assume the hardening and softening only depend on the volumetric plastic deformation \(\epsilon _V^P = {\mathrm {log}}\left( {{\mathrm {det}}\left( {{\bf{F}}^P} \right)} \right)\). We follow the derivation from Ortiz and Pandolfi35 and use the following hardening law:

$$p_0 = K{\kern 1pt} {\mathrm{sinh}}\left( {\xi {\kern 1pt} {\mathrm{max}}\left( { - \epsilon _V^P,0} \right)} \right),$$ (4)

where ξ is the hardening factor and K is the material bulk modulus. When the plastic deformation is compressive \(\left( {\dot \epsilon _V^P < 0} \right)\), p 0 will increase, causing the yield surface to grow in size. Snow will consequently receive more elastic responses resisting compression. When the plastic volume is increased \(\left( {\dot \epsilon _V^P > 0} \right)\), the yield surface shrinks which allows the snow to fracture in tension. This hardening law is represented in Fig. 1b (in black).

Classical hardening/softening laws such as the one described above for the dense snow slab (Eq. 4) fail in reproducing the collapse of porous cohesive materials under compression. This is shown in Fig. 1c, d (black lines) in which p significantly increases after reaching the yield surface and q slightly decreases before increasing. Hence, for the porous weak layer, we propose a modified softening law that describes cohesion and volume loss under compressive stresses. This new softening law involves looking at the volumetric plastic strain rate \(\dot \epsilon _V^P\). We introduce the anticrack plastic strain η which is related to \(\epsilon _V^P\) as follows:

$$\dot \eta = \left( {\begin{array}{*{20}{l}} {\alpha \left| {\dot \epsilon _V^P} \right|,} \hfill & {{\mathrm{if}}{\kern 1pt} t \le t_c} \hfill \\ {\dot \epsilon _V^P,} \hfill & {{\mathrm{if}}{\kern 1pt} t > t_c} \hfill \end{array}} \right.$$ (5)

where α is a softening factor which controls the energy dissipated during fracture and t c is the time corresponding to complete softening, i.e., \(\epsilon _V^P = 0\) and p 0 = 0 (state (2*) in Fig. 1). Our new softening law for the weak layer is obtained by replacing \(\epsilon _V^P\) by η in Eq. (4) (the discretization is shown in the Methods section). Hence, when stresses in the weak layer reach the yield surface, the introduction of the norm of \(\dot \epsilon _V^P\) in Eq. (5) will lead to softening (through a decrease in p 0 ) even under compression for which \(\dot \epsilon _V^P < 0\). The yield surface thus shrinks until it corresponds to a point at the origin of the p–q space. In addition, cohesion is removed by setting β = 0 when \(\epsilon _V^P = 0\) which ensures continuity. After reaching this point, the yield surface is free to expand according to the classical hardening law (Eq. 4), leading to volume reduction (collapse) due to the weight of the slab (blue arrows in Fig. 1b) and then a purely frictional/compaction behavior. Our softening rule reproduces bond breaking in the weak layer and subsequent grain rearrangement leading to volumetric collapse due to the compressive weight of the slab36. In contrast to classical hardening laws, our new formulation induces a strain-softening behavior even under macroscopic uniaxial compression, as shown in Fig. 1c, d. The observed mechanical behavior is very similar to that reported in discrete element simulations of porous cohesive granular materials37 and follows the following sequence of mechanical regimes: elastic regime, failure, drop in pressure and shear stress (strain softening), plastic consolidation corresponding to the volumetric collapse and, finally, dense packing regime corresponding to the jamming transition. This typical post-peak behavior was also observed in laboratory experiments of snow failure38,39 as well as during the propagation of compaction bands in confined compression of snow2. Physically, this behavior is related to the fact that even under a macroscopic compressive loading mode, the solid matrix of porous cohesive materials is mostly under tension (bending) and shear37. The behavior of the weak layer during a shear test simulation is shown in Supplementary Note 2.

To remove mesh dependency induced by softening, we follow the suggestion of Mahajan et al.20 and Sulsky and Peterson40 to regularize the jump in displacement. It is performed by dissipating the same amount of energy for different mesh resolutions by setting the softening factor α in Eq. (5) proportional to the mesh size dx. The influence of the mesh resolution on the volumetric plastic deformation during weak layer collapse is shown in Supplementary Note 3. For more detail about the Material Point Method and the implementation of the plastic model (plastic flow rule and return mapping), please see the Methods section.

Let us summarize here the different model parameters and their physical meaning: p 0 is the consolidation pressure and represents the compressive strength of the material, β is the ratio between tensile and compressive strength and represents cohesion, M is the slope of the Critical State Line (CSL) and characterizes the friction of the material, K is the bulk elastic modulus, ξ is the hardening coefficient and characterizes the brittleness of the material (a large ξ makes snow more brittle) and α is the softening factor which controls the fracture energy of the weak layer.

Field experiments

We report anticrack propagation in Propagation Saw Test (PST) experiments14. A PST consists in creating an artificial crack of increasing size by cutting within the weak layer with a saw until crack propagation. Depending on snowpack properties, the crack can either propagate until the end of the column (“END” case) or induce a fracture in the slab thus arresting the propagation (“SF” case). Black markers are inserted in the snowpack in order to derive the displacements using particle tracking velocimetry (PTV) and a high speed video camera. Two experiments were performed on flat terrain (ψ = 0°) and one on a typical avalanche slope (ψ = 37°). The density of the slab ranged from 159 to 279 kg m−3, slab thickness ranged from 26 to 75 cm and weak layer thickness ranged from 1 to 15 cm. For more detail about the experimental set-up, snowpack properties and data analysis, please refer to the Methods section.

Figure 2 shows the vertical displacements u y of the markers during the experiments as well as the displacement field at different key instants of the experiments. For each experiment, the crack in the weak layer induces slab bending leading to relatively small displacements. After reaching the critical crack length, the vertical displacement increased significantly due to dynamic anticrack propagation inducing the progressive collapse of the weak layer.

Fig. 2 Comparison between experimental and simulated results. Experimental (top) and numerical (bottom) results for experiment number 1 (a), experiment number 2 (b), and experiment number 3 (c). The displacement field is shown on the left for different key instants in each case: during anticrack propagation (t = 6.65 s) in a; during frictional sliding (t = 3.7 s) in b, and after slab fracture (t = 2.5 s) in c. On the right, the time evolution of the average vertical displacement of vertical rows of markers is shown. The color represent the average horizontal position of each vertical row of markers. ACP anticrack propagation. The red color in the weak layer represents plasticity Full size image

The first experiment on the flat (experiment number 1) highlights the potential of remote avalanche triggering from low-angle terrain. All markers show significant vertical displacements (between 2.5 and 8 mm) and the fracture in the weak layer propagated until the end of the beam (END). The second experiment (experiment number 2) made on a typical avalanche slope is also a case of full propagation (END) with significant collapse of the weak layer (up to 1 cm). In this case, crack propagation is followed by the sliding of the slab since slope angle is larger than the friction angle of snow (~30°, van Herwijnen et al.41). Sliding induces the progressive erosion of the weak layer and thus further vertical displacement. Full propagation in the weak layer is typical for deep and dense slab layers16. The third experiment (experiment number 3) on the flat is a case of partial propagation in the weak layer with slab fracture (SF). In this case, markers located on the right side of the fracture did not move. This is a typical outcome for low density and shallow slab layers15,16,42. Movies of the three experiments, including displacement fields are provided in the Supplement (Supplementary Movies 1–3).

PST simulations

The hybrid Eulerian–Lagrangian Material Point Method was used to solve the set of partial differential equations of the system, given the same characteristics and boundary conditions as in the experimental PSTs. We discuss in detail the choice of snowpack mechanical properties in the Methods section.

As shown in Fig. 2 and in the Supplementary Movies, our model accurately reproduces all the features observed in the experiments. More specifically, anticrack propagation on flat terrain i.e. without external driving shear forces, is very well captured (Fig. 2a and Supplementary Movie 1). A measured critical crack length a c = 39 cm was well reproduced by the simulation. The collapse wave speed c was computed from the time-delay between the onset of movement between markers15. It was around 35 m s−1 in both the experiment and the simulation. This speed is significantly lower than the speed of elastic waves \(c^e = \sqrt {E{\mathrm{/}}\rho }\) in the slab which is around 200 m s−1. Once the crack has propagated through the full system length, the system is at rest. Figure 2b and Supplementary Movie 2 show the results on a typical avalanche slope of 37°. Anticrack propagation features are very similar as on the flat but the propagation speed is lower (c = 23 m s−1) due to a lower slab elastic modulus and density and a larger weak layer strength16 than in experiment number 1. The bending phase, critical crack length (a c = 32 cm), anticrack propagation as well as the frictional sliding of the slab are very well reproduced by our model. Crack branching resulting from the interplay between weak layer and slab fracture is also well reproduced, as shown on Fig. 2c and Supplementary Movie 3. In this case, anticrack propagation in the weak layer was arrested 10 cm after reaching a critical length (a c = 26.5 cm) as the tensile stress in the slab induced by slab deformation exceeded the tensile strength due to its thin and weak character (low density slab). In contrast, in the two previous experiments, the tensile stress in the slab remained lower than the strength thus leading to full propagation. Nevertheless, we note that the bending deformation pre-propagation was underestimated by our model for experiment number 3. This suggests that inelastic (probably rate-dependent) deformation is induced by the very loose character of the slab in this experiment (ρ 3 = 159 kg m−3). In addition, we observe small oscillations in our displacements because our simulations are performed without damping.

Note that for all simulations, the anticrack velocity was found almost equal to the collapse speed obtained from the vertical displacement of the slab. However, we observed that the anticrack tip is always located slightly ahead of the collapse wave front.

Slope-scale simulations

Two- and three-dimensional slope-scale simulations of remote avalanche triggering were performed (Supplementary Movies 4–7). In both 2D and 3D slope simulations, the average crack propagation speed was around 60 m s−1 and the crown fracture was almost perpendicular to the bed surface as reported by Perla43 and McClung and Schweizer44. Furthermore, the slab fracture at the crown of the avalanche (upslope section of the fracture line) started branching from the bottom of the slab at the interface with the weak layer (Supplementary Movie 5), in contrast to the PST simulation and experiment number 3 in which it started branching from the top. In 3D, the simulated release zone (Fig. 3, Supplementary Movie 7) has commonly observed characteristics43: an arc crown line as well as jagged flanks (side sections of the fracture line) and staunchwall (bottom section of the fracture line). Crown fracture occurs in tension while flank and staunchwall fractures occur in shear. Finally, the cross-slope propagation was approximately twice slower than up-slope propagation.