the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
Micromechanical modeling of snow failure
Grégoire Bobillier
Bastian Bergfeld
Achille Capelli
Jürg Dual
Johan Gaume
Alec van Herwijnen
Jürg Schweizer
Drysnow slab avalanches start with the formation of a local failure in a highly porous weak layer underlying a cohesive snow slab. If followed by rapid crack propagation within the weak layer and finally a tensile fracture through the slab, a slab avalanche releases. While the basic concepts of avalanche release are relatively well understood, performing fracture experiments in the laboratory or in the field can be difficult due to the fragile nature of weak snow layers. Numerical simulations are a valuable tool for the study of micromechanical processes that lead to failure in snow. We used a threedimensional discrete element method (3D DEM) to simulate and analyze failure processes in snow. Cohesive and cohesionless ballistic deposition allowed us to reproduce porous weak layers and dense cohesive snow slabs, respectively. To analyze the micromechanical behavior at the scale of the snowpack (∼1 m), the particle size was chosen as a compromise between low computational costs and detailed representation of important micromechanical processes. The 3DDEM snow model allowed reproduction of the macroscopic behavior observed during compression and mixedmode loading of drysnow slab and the weak snow layer. To be able to reproduce the range of snow behavior (elastic modulus, strength), relations between DEM particle and contact parameters and macroscopic behavior were established. Numerical loadcontrolled failure experiments were performed on small samples and compared to results from loadcontrolled laboratory tests. Overall, our results show that the discrete element method allows us to realistically simulate snow failure processes. Furthermore, the presented snow model seems appropriate for comprehensively studying how the mechanical properties of the slab and weak layer influence crack propagation preceding avalanche release.
 Article
(1821 KB)  Fulltext XML

Supplement
(485 KB)  BibTeX
 EndNote
Drysnow slab avalanches require initiation and propagation of a crack in a weak snow layer buried below cohesive slab layers. Crack propagation occurs if the initial zone of damage in the weak layer is larger than the socalled critical crack size. Weak layer fracture during crack propagation is generally accompanied by the structural collapse of the weak layer due to the high porosity of snow (van Herwijnen et al., 2010). If the crack propagates across a steep slope, a slab avalanche may release (McClung, 1979; Schweizer et al., 2003). Our understanding of crack propagation was greatly improved by the introduction of the propagation saw test (PST; Gauthier and Jamieson, 2006; Sigrist and Schweizer, 2007; van Herwijnen and Jamieson, 2005). The PST involves isolating a snow column and initiating a crack by sawing in a predefined weak layer until the critical crack length is reached and selfpropagation starts. The PST allows analysis of the onset and dynamics of crack propagation and derivation of important mechanical properties using particle tracking velocimetry (van Herwijnen et al., 2016). The essential mechanical properties related to the onset of crack propagation are slab elasticity, slab load and tensile strength, as well as the weak layer strength and specific fracture energy (e.g., Reuter and Schweizer, 2018). However, no theoretical framework exists that describes how these mechanical properties and possibly other ones such as the weak layer failure envelope, weak layer elasticity or microstructure relate to the dynamics of crack propagation at the slope scale. Whereas field experiments are difficult to perform at this scale, numerical simulations may provide insight into the drivers of propagation dynamics.
Johnson and Hopkins (2005) were the first to apply the discrete element method (DEM) to model snow deformation. They simulated creep settlement of snow samples, which consisted of 1000 randomly oriented cylinders of random length with hemispherical ends. More recently, the DEM was used to model the mechanical behavior based on the complete 3D microstructure of snow (Hagenmuller et al., 2015). Gaume et al. (2015) developed a discrete element model to simulate crack propagation and subsequently derived a new analytical expression for the critical crack length (Gaume et al., 2017b). Their approach allows the generation of highly porous samples and was used to perform 2D simulations of the PST in agreement with field experiments. However, the oversimplified shape (triangular structure) and the 2D character of the weak layer employed by Gaume et al. (2015) prevented a detailed analysis of the internal stresses during crack propagation. On the other hand, microstructurebased DEM models adequately reproduce the mechanical behavior (Mede et al., 2018). However, the computational costs of these complex 3D models are too high to generate samples large enough to investigate the dynamics of crack propagation at the slope scale on a standard personal computer (Intel Core i7, eight processors, 3.4 GHz, RAM 16Gb).
Our aim is to develop a 3DDEM snow model that adequately takes into account snow microstructure, but is not too costly in terms of computational power so that simulations at the slope scale become feasible on a highperformance computer (HPC). To relate DEM parameters to macroscopic snow behavior, we will validate the model by simulating basic load cases. Finally, we numerically simulate mixedmode loading experiments and compare results to those obtained during laboratory experiments.
2.1 Formulation of the model
2.1.1 Discrete element method
To simulate the failure behavior of layered snow samples, we used the discrete element method (DEM). DEM, first introduced by Cundall and Strack (1979), is a numerical tool, commonly employed to study granularlike assemblies composed of a large number of discrete interacting particles. We used the PFC3D (v5) software developed by Itasca Consulting Group (http://www.itascacg.com, last access: 20 December 2019).
2.1.2 Contact model
We used the parallelbond contact model (PBM) introduced by Potyondy and Cundall (2004). The PBM provides the mechanical behavior of a finitesized piece of cementlike material that connects two particles. The PBM component acts in parallel with a classical linear contact model and establishes an elastic interaction between the particles. The mechanical parameters include the contact elastic modulus E_{u}, Poisson's ratio ν_{u}=0.3, restitution coefficient e_{u}=0.1 and friction coefficient ${\mathit{\mu}}_{\mathrm{u}}=\mathrm{0.2}.$ If particles are bonded, the bond part will act in parallel to the contact part. The bonded part is described by the bond elastic modulus E_{b}, the bond Poisson's ratio ν_{b}=0.3 and the bond strength, shear and tensile strength σ_{s} and σ_{t}. To reduce the number of variables we assume E_{u}=E_{b}≜E_{particle} and ${\mathit{\sigma}}_{\mathrm{s}}={\mathit{\sigma}}_{\mathrm{t}}\mathit{\triangleq}{\mathit{\sigma}}_{\mathrm{bond}}^{\mathrm{th}}$. Figure 1 illustrates the PFC parallelbond model (PBM) with the mechanical parameters for the bonded and unbonded states. Four different bond behaviors (tension–compression, shear, bending and torsion) are shown in Fig. 2. The maximum bond shear and tensile stresses are calculated via beam theory depending on the normal and shear forces in the bond, M_{b,1} the bending moment, M_{b,2} the twisting moment, r_{b} the bond radius, A the area of the bond, I the moment of inertia of the bond cross section and J the polar moment of inertia. More details on the PBM can be found in previous studies (Gaume et al., 2015, 2017a, b). Once a bond breaks, only particle frictional contact occurs, and no new bonds are created (i.e., no sintering occurs). This assumption is motivated by the fact that the strain rate is large and the timescale is seconds during the postfailure phase.
2.1.3 System generation
The simulated threedimensional system consisted of a rigid basal layer (Fig. 3, blue particles), the layer studied (weak layer or slab layer, green particles), and an “actuator” layer used to apply the load (red particles). The basal layer is composed of a single layer of particles with a radius of r=5 mm. The weak layer was created by cohesive ballistic deposition (Löwe et al., 2007) to reproduce the porous and anisotropic structure of natural weak layers. Doing so, we obtained a porosity of 80 % for a particle radius of r=2.5 mm. The layer thickness (3 cm) can be modified by homothetic transformation while keeping the same mechanical behavior. A short weak layer scaling study is provided in the Supplement.
We used cohesionless ballistic deposition to generate dense layers (Kadau and Herrmann, 2011) as typically found in snow slab layers. For these layers we used a particle radius of $r=\mathrm{11}\pm \mathrm{1}$ mm (uniform distribution). The radius variation was introduced to prevent close packing, resulting in a porosity of 45 %. Layer density (ρ) was adjusted by changing the particle density. The size of the particles is not intended to represent the real snow grains. The particle size was chosen as a tradeoff between an acceptable computation time (minutes to days) and avoiding particle size effects in the numerical experiments. At the defined particle scale (larger than the snow grains) the ice properties (e.g., strength, elastic modulus, Poisson's ratio) cannot be used directly. Therefore, the particle scale can be considered a mesoscale between the macroscopic scale (sample scale) and the microscale (individual snow grains). Hence, we adjusted the particle density to represent the macroscopic snow densities in accordance with the macroscopic sample porosity.
To characterize the mechanical behavior of these two types of snowlike layers (weak layer or slab layer), unconfined loadcontrolled tests were performed and compared to experimental results. We also performed confined compression tests, but we found no difference in behavior compared to the unconfined tests due to the porosity of 80 % (shown in the Supplement). To simulate the tests, we added an “actuator” layer generated by cohesionless ballistic deposition, composed of particles with a radius r=10 mm on top of the studied layer (Fig. 3, red particles). This layer is defined as a rigid clump with initially low density. A clump is a rigid collection of n rigid particles that form one DEM element. The volume is defined by the particle positions and radius. The mass properties are defined by the clump density and clump volume. Clumps can translate and rotate but cannot deform. Clump motion obeys the equations of motion induced by the definition of mass properties, loading conditions and velocity conditions.
The samples were generated in a box; the box walls were then removed to create unconfined test conditions. To avoid a packing effect at the sidewalls, samples were generated 10 particle radii larger and cutout before the simulation. In order to model macroscopic mechanical behavior of the studied layers, we tuned the particle elastic moduli and the bond strength. A large range of particle elastic modulus and bond strengths were tested to characterize the relation between particle parameters and macroscopic behavior. In some materials, strength and elastic modulus are related, while in other materials these properties are independent. For snow, it remains unclear whether the two properties are related. Our goal was to independently control both parameters in order to have a precise control on the macroscopic elastic modulus and macroscopic strength of the snow.
2.1.4 Loadcontrolled test
Loadcontrolled simulations were performed by linearly increasing the actuator layer density. To avoid a sample size effect (see below), 30 cm×30 cm samples were generated. Our DEM model does not take into account viscous effects or sintering of snow; therefore the results do not depend on the loading rate. We chose a high loading rate of 20 kPa s^{−1} simply to reduce the simulation time but we verified that the loading rate did not affect the results (not shown). By changing the angle of gravity (ψ in Fig. 3), mixedmode loading was simulated. Failure was defined as the point of maximum shear or normal stress during the two first phases (linear elastic and softening phases). Table 1 summarizes the particle mechanical properties used for simulations and their corresponding macroscopic values.
2.1.5 Time step
The length of the time step was determined as a function of the particle properties according to
where ρ and r are the smallest particle density and radius, respectively, and E is the largest bond or particle elastic modulus. Choosing the time step in this manner ensures the stability of the DEM model (Gaume et al., 2015).
2.1.6 Stress and strain
The average stress and strain were calculated at the interface between the rigid base layer and the studied layer (Fig. 3, violet arrow). Normal stress ${\stackrel{\mathrm{\u203e}}{\mathit{\sigma}}}_{z}$ was computed as ${\stackrel{\mathrm{\u203e}}{\mathit{\sigma}}}_{z}={F}_{z}/A$ and shear stress as ${\stackrel{\mathrm{\u203e}}{\mathit{\sigma}}}_{x}={F}_{x}/A$. Here F_{x} and F_{z} are the sum of the contact forces acting on the basal layer in the tangential and normal directions, respectively, and A is the total area of the basal layer over which the stresses were determined. We define the engineering strains as normal strain, ${\mathit{\epsilon}}_{z}=\frac{{u}_{z}}{D}$, and shear strain, ${\mathit{\gamma}}_{x}=\frac{{u}_{x}}{D}$, with the displacement of the actuator u in the z and x directions and the thickness D of the studied layer. The macroscopic strength (σ^{th}) was defined as the maximum stress before catastrophic failure. The macroscopic elastic modulus (E) was defined on the normal stress–strain curve as the derivative of the stress between 0 % and 70 % of the stress peak.
2.1.7 Fabric tensor
If the particle arrangement during layer creation is not isotropic, the mechanical quantities of the layer show directional dependency. For any heterogeneous, anisotropic material (e.g., bones, concrete, snow), the fabric tensor characterizes the geometric arrangement of the porous material microstructure. The fabric tensor, referred to here as the contact tensor C, is the volume average of the tensor product of the contact unit normal vectors $\stackrel{\mathrm{\u203e}}{n}$. The 2ndorder contact tensor coefficients are defined in KenIchi (1984) as
where N is the total number of particle contacts, and ${n}_{i}^{\mathit{\alpha}}$ represents the normalized projections of the contact with respect to the x_{i} Cartesian coordinate (Shertzer et al., 2011). The contact tensor C was used to estimate the physical properties of the simulated sample.
2.1.8 Laboratory experiments
For model validation, we used data of cold laboratory experiments obtained with a loading apparatus described in Capelli et al. (2018). They performed loadcontrolled failure experiments on artificially created, layered snow samples, consisting of a weak layer of depth hoar crystals between harder layers of finegrained snow. The load applied to the samples was increased linearly until the sample failed. For more information on the experiments, see Capelli et al. (2018). We selected three experiments (Table 2) for validating the numerical simulations. For the validation we focused on the normal strain, since for the experimental shear strain data (measure of the horizontal displacement) the signaltonoise ratio was too low. Furthermore, due to the method used to load the snow samples, data from the force sensor after failure contained experimental artifacts. To select the model parameters E_{particle} and ${\mathit{\sigma}}_{\mathrm{bond}}^{\mathrm{th}}$, we used the elastic modulus computed as the derivative of the normal stress–strain curve and the strength values from the experiments (Table 2), as well as the relations for strength and modulus derived below. Digital image analysis of the experiments had revealed that the deformation was concentrated in the weak layer (Capelli et al., 2018). We therefore simulated the weak layer with a rigid actuator layer on the top.
This section first presents the structural properties of the two generated layers. The two generated layers were analyzed based on an unconfined compression test. Then, the link between macroscopic behavior and particle properties is described. Finally, the model setup for the weak layer is validated by comparing numerical mixedmode loading simulations to experimental data.
3.1 Structural properties of generated samples
For the sample used to generate the slab, the coefficients of the contact tensor C were (Eq. 2)
This shows that the slab samples are nearly isotropic, which is in line with results reported for snow types typically found in snow slab layers (Gerling et al., 2017; Srivastava et al., 2016).
For the weak layer samples, five in total, the contact tensor was
It shows transverse isotropic symmetry that is again in line with data from snow samples representative for weak layers (i.e., layers of depth hoar, surface hoar or facets), which also show transverse isotropy (Gerling et al., 2017; Shertzer, 2011; Shertzer et al., 2011; Srivastava et al., 2016).
3.2 Characterization of macroscopic properties
3.2.1 Slab layer
To establish a relationship between the macroscopic elastic modulus and the particle elastic modulus, we performed 100 simulations (with 10 different values of E_{particle} and 10 different values of ${\mathit{\sigma}}_{\mathrm{bond}}^{\mathrm{th}})$ in pure compression to relate macroscopic and particle parameters. The macroscopic elastic modulus increased linearly with E_{particle}:
with the coefficients ${\mathit{\beta}}_{\mathrm{0}}=\mathrm{1.5}\times {\mathrm{10}}^{\mathrm{5}}$ Pa and β_{1}=0.526 (dashed line in Fig. 4a).
The macroscopic strength also increased linearly with bond strength:
with the coefficients ${\mathit{\gamma}}_{\mathrm{0}}=\mathrm{318}$ Pa and γ_{1}=0.982 (dashed line in Fig. 4b; R^{2}=0.999).
3.2.2 Weak layer
For the weak layer we performed 81 simulations (with nine different values of E_{particle} and nine different values of ${\mathit{\sigma}}_{\mathrm{bond}}^{\mathrm{th}})$ in pure compression to relate macroscopic and particle parameters. The macroscopic elastic modulus increased linearly with E_{particle}:
with coefficients ${\mathit{\beta}}_{\mathrm{0}}=\mathrm{7.3}\times {\mathrm{10}}^{\mathrm{4}}$ Pa and β_{1}=0.014 (Fig. 5a; R^{2}=0.985).
The macroscopic strength also increased linearly with bond strength:
with coefficients γ_{0}=76.7 Pa and γ_{1}=0.016 (Fig. 5b; R^{2}=0.998).
Hence, based on Eqs. (7) and (8), any values of the macroscopic quantities ${\mathit{\sigma}}_{\mathrm{wl}\phantom{\rule{0.25em}{0ex}}\mathrm{macro}}^{\mathrm{th}}$ and E_{wl macro} can be obtained by tuning E_{particle} and ${\mathit{\sigma}}_{\mathrm{bond}}^{\mathrm{th}}$.
3.3 Mechanical behavior of layers
3.3.1 Slab layer
To investigate the mechanical behavior of the slab layer, we performed loadcontrolled tension tests. Two phases can be distinguished: linear elastic deformation followed by sample fracture. During the linear elastic deformation, no bond damage appears, and the stress linearly increases up to ε=0.0025 (Fig. 6). At failure, the stress dropped rapidly and bond damage drastically increased with increasing strain.
3.3.2 Weak layer
The largedeformation, unconfined loadcontrolled compression tests of weak layer samples revealed four different phases (grey dashed–dotted lines in Fig. 7). First, there was a linear elastic phase without bond breaking (a.1); nonlinearity appears right before the stress peak induced by some damage prior to failure, in good agreement with the quasibrittle behavior of weak snow layers (Fig. 8). When the macroscopic strength was reached, the normal stress dropped sharply during the softening phase as bond damage increased drastically (a.2). During the brittle crushing phase, the sample density as well as the proportion of broken bonds (P_{broken bond}) steadily increased (a.3). Finally, the densification phase (a.4) was reached when the stress prominently increased again as the particles were closely packed.
By varying the particle modulus E_{particle} and the bond strength ${\mathit{\sigma}}_{\mathrm{bond}}^{\mathrm{th}}$, the micromechanical behavior in terms of bond breaking and acceleration (a) of the actuator layer was also investigated more closely up to the start of the brittle crushing phase (Fig. 9). Before reaching the macroscopic strength, the normal stress increased linearly with increasing strain while the number of broken bonds and the acceleration were low. The strain at failure depended on both E_{particle} and ${\mathit{\sigma}}_{\mathrm{bond}}^{\mathrm{th}}$. During the softening the stress sharply dropped while both the number of broken bonds and the acceleration increased. Both E_{particle} and ${\mathit{\sigma}}_{\mathrm{bond}}^{\mathrm{th}}$ controlled the amount of stress drop as well as the rate of increase in P_{broken bond} and a. During the brittle crushing phase, both strength and acceleration did not change while P_{broken bond} increased, independent of the values of E_{particle} and ${\mathit{\sigma}}_{\mathrm{bond}}^{\mathrm{th}}$.
The stress at the end of the softening phase was characterized by the softening ratio $R=\frac{{\widehat{\mathit{\sigma}}}_{\mathrm{residual}}}{{\mathit{\sigma}}_{\mathrm{wl}\phantom{\rule{0.25em}{0ex}}\mathrm{macro}}^{\mathrm{th}}}$, with ${\mathit{\sigma}}_{\mathrm{wl}\phantom{\rule{0.25em}{0ex}}\mathrm{macro}}^{\mathrm{th}}$ the macroscopic strength and ${\widehat{\mathit{\sigma}}}_{\mathrm{residual}}$ the mean residual stress during the brittle crushing phase. The test with the highest softening ratio (Fig. 9 solid light blue line: R=0.45) showed the lowest damage and the lowest acceleration. In contrast, the lowest softening ratio (Fig. 9 dark blue dashed line: R=0.21) corresponded to the largest proportion of broken bonds and the largest acceleration. Concerning the two other tests, they exhibited the same residual stress but different softening ratios. We observed that the softening ratio followed a nonlinear relation with E_{particle} and ${\mathit{\sigma}}_{\mathrm{bond}}^{\mathrm{th}}$.
Similar to the behavior under compression, the mechanical response in shear exhibited different phases: an elastic phase, softening and simultaneously normal brittle crushing and shear displacement, and finally shear displacement only (grey dashed–dotted lines in Fig. 10). Also, the damage dynamics were similar to in pure compression (Fig. 10b). No critical bond breaking was observed during the linear elastic phase followed by catastrophic damage after failure. Subsequently, the damage further increased during the brittle crushing. The normal strain increased during the brittle crushing phase and did not change thereafter. The normal deformation was closely related to the proportion of broken bonds, similar to behavior in the pure compression. Shear and normal accelerations reached their maximum at the end of the softening phase (Fig. 10c) as observed in pure compression (Fig. 7). During the brittle crushing phase, the normal acceleration decreased due to the creation of new contacts that decelerate the actuator layer. The tangential acceleration did not change much during the shear displacement phase.
3.3.3 Weak layer failure envelope
Unconfined loadcontrolled tests with nine loading angles were performed to create the failure envelope. Figure 11 compiles the values of macroscopic strength for different loading angles resulting in a failure envelope including tension (negative normal stress), pure shear, pure compression and mixedmode loading states. To investigate the influence of sample size, we performed a sensitivity analysis by varying the sample size from 0.1 m×0.1 m to 1 m×0.6 m and the random deposition (generation of different ball positions for the ballistic deposition). Apart from the smallest sample, all samples had very similar failure envelopes, which were fitted with a 2ndorder polynomial with coefficients ${\mathit{\beta}}_{\mathrm{0}}^{\mathrm{FE}}=\mathrm{7.66}\times {\mathrm{10}}^{\mathrm{2}}\phantom{\rule{0.125em}{0ex}}\mathrm{Pa}$, ${\mathit{\beta}}_{\mathrm{1}}^{\mathrm{FE}}=\mathrm{0.25}$ and ${\mathit{\beta}}_{\mathrm{2}}^{\mathrm{FE}}=\mathrm{1.97}\times {\mathrm{10}}^{\mathrm{4}}$ (dashed–dotted black line in Fig. 11, R^{2}=0.97):
For a sample length of 0.3 m or larger, no effect of sample size on the failure envelope was observed. The sample heterogeneity induced by different types of random deposition did not influence the failure envelope either. Given the expression for the macroscopic strength (Eq. 9), where σ^{th} represents the normal strength and τ^{th} the shear strength, the failure envelope is directly related to ${\mathit{\sigma}}_{\mathrm{bond}}^{\mathrm{th}}$.
As the macroscopic strength ${\mathit{\sigma}}_{\mathrm{wl}}^{\mathrm{th}}$ is related to ${\mathit{\sigma}}_{\mathrm{bond}}^{\mathrm{th}}$ (Eq. 8), the failure envelope can be scaled by using the scaling factor ($\frac{{\mathit{\sigma}}_{\mathrm{wl}}^{\mathrm{th}}}{{\mathit{\sigma}}_{\mathrm{wl}\phantom{\rule{0.25em}{0ex}}\mathrm{ref}}^{\mathrm{th}}})$:
with ${\mathit{\sigma}}_{\mathrm{wl}\phantom{\rule{0.25em}{0ex}}\mathrm{ref}}^{\mathrm{th}}=\mathrm{2650}\phantom{\rule{0.125em}{0ex}}\mathrm{Pa}$, which corresponds to the maximum strength in pure compression (Fig. 11). Equation (10) allows the derivation of the failure envelope for any value of the bond strength ${\mathit{\sigma}}_{\mathrm{bond}}^{\mathrm{th}}$ (green dashed–dotted lines in Fig. 12).
3.4 Comparison with experimental data
To validate the behavior of our simulated weak layer samples, we used data from laboratory experiments performed by Capelli et al. (2018) (Table 2). For each of the three experiments with different loading angles, the simulated total stress (${\mathit{\sigma}}_{\mathrm{tot}}=\sqrt{{\mathit{\sigma}}^{\mathrm{2}}+{\mathit{\tau}}^{\mathrm{2}}})$ as a function of normal strain (ε) was in good agreement with the experimental results (Fig. 13).
We used 3D discrete element modeling to study the mechanical behavior of simplified snow samples generated by different ballistic deposition techniques. Cohesive ballistic deposition produced transversally isotropic weak layers with high porosity (80 %). Cohesionless ballistic deposition produced isotropic slab layers of lower porosity (45 %), in general agreement with key properties of natural snow samples (Shertzer, 2011). The DEM particles do not represent real snow grains to keep the computational costs reasonable (i.e., ∼10 min on a standard personal computer for a sample of 50 cm×50 cm in size, corresponding approximately to 26 500 particles). By varying the DEM particle parameters E_{particle} and ${\mathit{\sigma}}_{\mathrm{bond}}^{\mathrm{th}}$, the macroscopic properties can be modified to fit different types of snow.
First, tension tests were simulated to study the behavior of dense slab layers. The results evidenced an almost perfectly brittle behavior (Fig. 6) in good agreement with the tensile behavior reported by Hagenmuller et al. (2014) and by Sigrist (2006).
The mechanical behavior we observed for our weak layer samples, in particular the four phases (Fig. 7) during a loadcontrolled compression test, was very similar to that reported by Mede et al. (2018), who simulated snow behavior with microstructurebased snow samples. More generally, Gibson and Ashby (1997) also described these four distinct phases for elastic–brittle foam samples.
The unconfined loadcontrolled tests under mixedmode loading conditions showed shear behavior in good agreement with previously reported results (Mede et al., 2018; Mulak and Gaume, 2019; Reiweger et al., 2015).
The obtained failure envelopes were qualitatively in good agreement with the Mohr–Coulomb–cap (MCC) model proposed by Reiweger et al. (2015) and with the ellipsoid (cam clay) model proposed by Gaume et al. (2018) and Mede et al. (2018). The model qualitatively reproduced the snow failure envelopes found in other numerical studies (Mede et al., 2018; Mulak and Gaume, 2019). In our case, the failure envelope is directly linked to ${\mathit{\sigma}}_{\mathrm{bond}}^{\mathrm{th}}$, since any failure envelope can be expressed as a function of ${\mathit{\sigma}}_{\mathrm{bond}}^{\mathrm{th}}$. Weak layer failure behavior was not affected by the heterogeneity induced by different types of random ball deposition and by the sample size if the sample size was larger than 0.3 m×0.3 m. This size is typically found in field tests (PST, extended column test; Bair et al., 2014; Reuter et al., 2015; van Herwijnen et al., 2016) and laboratory experiments (Capelli et al., 2018).
Based on these purely numerical investigations, the particle and contact parameters were selected to reproduce the results of cold laboratory experiments with real snow samples (Fig. 13). The numerical results were qualitatively in good agreement for the three loading angles. However, the comparison to the experimental results is hindered by the lack of adequate experimental data. Due to vibrations in the actuator plate, the experimental shear strain data could not be used. Hence, there are no experimental data to validate the postfailure behavior. Still, as shown above, the postfailure behavior was in agreement with results of other numerical studies (e.g., Mede et al., 2018).
We showed that the onset of failure corresponded to a strong increase in the number of broken bonds and in actuator layer acceleration (Fig. 9). The maximum acceleration was reached towards the end of the softening phase. In fracture mechanics, the zone where softening occurs is generally referred to as the fracture process zone (FPZ) (Bazant and Planas, 1998). Hence, our findings suggest that slab acceleration may be used to accurately track the crack tip location in the weak layer during crack propagation experiments.
Introducing the softening ratio (R) showed that the stress decrease in softening only depends on particle modulus E_{particle} and bond strength ${\mathit{\sigma}}_{\mathrm{bond}}^{\mathrm{th}}$, which allows estimation of the maximum acceleration of the actuator layer and the damage dynamics. In the present formulation of our model, the softening ratio is fixed for a given pair of parameters (E_{particle} and ${\mathit{\sigma}}_{\mathrm{bond}}^{\mathrm{th}})$.
To limit the number of model parameters we made two assumptions: the contact and the bond elastic moduli are equal and the bond cohesive and tensile strengths are equal. The choice of weak layer creation technique (cohesive ballistic deposition) caused unique structural anisotropy that was reflected in the mechanical behavior and added a limitation to the postfailure behavior and the shape of the failure envelope. Investigating the influence of microstructure by modifying the porosity or the coordination number as the sticky hard sphere (Gaume et al., 2017a) and/or modifying the assumption on contact/bond elastic moduli would allow us to generate a larger range of stress decrease during the softening phase.
Furthermore, in the future, the influence of the ratio between the bond tensile strength and the bond cohesive strength, and/or the weak layer microstructure on the yield surface might be explored.
The developed simulation tool does not take into account snow sintering processes, as we limited the study to fast loading rates. In the context of a drysnow slab avalanche formation, this means that we can only study artificially induced cracks due to skiers or explosives. In the future, we plan to extend the work to larger systems with the objective of studying the micromechanics of the dynamics of crack propagation. Using the presented tool to model a PST already showed some promising preliminary results (Bobillier et al., 2018).
Understanding the failure behavior of slab and weak layer independently and characterizing the influence of the main parameters are prerequisites for studying the dynamics of crack propagation leading to the release of a drysnow slab avalanche.
We developed a mesoscale (between snow grain and slope scale) simulation tool based on 3D discrete element simulations to generate snow layers of varying properties and investigate micromechanical processes at play during snow failure. Two types of snow layers were generated using a ballistic deposition technique: (1) a uniform snow slab and (2) a porous transversally isotropic weak snow layer. These two types of snow layers are the two main components of drysnow slab avalanches. The layers were characterized by a linear relation between particle and contact parameters and macroscopic properties. By deliberately making the choice of not representing the real snow microstructure, the computational time decreases and allows the creation of relatively large systems.
We found an elastic–brittle mechanical behavior for slab layers in tension. The weak layer behavior under mixedmode loading included four distinct phases of deformation (elastic, softening, simultaneous brittle crushing and shear displacement, and finally shear displacement) as recently reported in the literature. The weak layer failure envelope, derived from a series of mixedmode loading simulations under different loading angles, was in good agreement with previous experimental and numerical results. The closedform failure envelope can be tuned by adjusting the bond strength parameter.
Analyzing weak layer features such as the proportion of broken bonds, normal acceleration and softening ratio showed some of the limitations induced by our assumptions on particle parameters and the uniqueness of the microstructure generation. Still, the validation results suggest that the presented simulation tool can reproduce the main behavior of weak layers under mixedmode loading conditions – even though we strongly simplified the microstructure to limit the computational costs.
In the future, we intend to increase the system size and simulate a propagation saw test and explore the dynamics of crack propagation that eventually lead to drysnow slab avalanche release.
All relevant data are available at https://doi.org/10.16904/envidat.122 (Bobillier et al., 2019).
The supplement related to this article is available online at: https://doi.org/10.5194/tc14392020supplement.
JS and AH initiated this study. GB developed the model and made all the simulations with help from JG. AC provided the laboratory data for validation. GB prepared the paper with contributions from all coauthors.
The authors declare that they have no conflict of interest.
We thank the anonymous reviewer, the reviewer Chris Borstad and the editor Evgeny A. Podolskiy for their constructive reviews that helped us to improve our paper.
Grégoire Bobillier has been supported by a grant from the Swiss National Science Foundation (200021_16942).
This paper was edited by Evgeny A. Podolskiy and reviewed by Chris Borstad and one anonymous referee.
Bair, E. H., Simenhois, R., van Herwijnen, A., and Birkeland, K.: The influence of edge effects on crack propagation in snow stability tests, The Cryosphere, 8, 1407–1418, https://doi.org/10.5194/tc814072014, 2014.
Bazant, Z. P. and Planas, J.: Fracture and size effect in concrete and other quasibrittle materials, CRC Press, Boca Raton, USA, 616 pp., 1998.
Bobillier, G., Gaume, J., van Herwijnen, A., Dual, J., and Schweizer, J.: Modeling the propagation saw test with discrete elements, Proceedings ISSW 2018, International Snow Science Workshop, Innsbruck, Austria, 7–12 October 2018, 976–980, 2018.
Bobillier, G., Bergfeld, B., Capelli, A., Dual, J., Gaume, J., van Herwijnen, A., and Schweizer, J.: Modeling snow failure with DEM, EnviDat, https://doi.org/10.16904/envidat.122, 2019.
Capelli, A., Reiweger, I., and Schweizer, J.: Acoustic emissions signatures prior to snow failure, J. Glaciol., 64, 543–554, https://doi.org/10.1017/jog.2018.43, 2018.
Cundall, P. A. and Strack, O. D. L.: A discrete numerical model for granular assemblies, Geotechnique, 29, 47–65, https://doi.org/10.1680/geot.1979.29.1.47, 1979.
Gaume, J., van Herwijnen, A., Chambon, G., Birkeland, K. W., and Schweizer, J.: Modeling of crack propagation in weak snowpack layers using the discrete element method, The Cryosphere, 9, 1915–1932, https://doi.org/10.5194/tc919152015, 2015.
Gaume, J., Löwe, H., Tan, S., and Tsang, L.: Scaling laws for the mechanics of loose and cohesive granular materials based on Baxter's sticky hard spheres, Phys. Rev. E, 96, 032914, https://doi.org/10.1103/PhysRevE.96.032914, 2017a.
Gaume, J., van Herwijnen, A., Chambon, G., Wever, N., and Schweizer, J.: Snow fracture in relation to slab avalanche release: critical state for the onset of crack propagation, The Cryosphere, 11, 217–228, https://doi.org/10.5194/tc112172017, 2017b.
Gaume, J., Gast, T., Teran, J., van Herwijnen, A., and Jiang, C.: Dynamic anticrack propagation in snow, Nat. Commun., 9, 3047, https://doi.org/10.1038/s4146701805181w, 2018.
Gauthier, D. and Jamieson, J. B.: Towards a field test for fracture propagation propensity in weak snowpack layers, J. Glaciol., 52, 164–168, https://doi.org/10.3189/172756506781828962, 2006.
Gerling, B., Löwe, H., and van Herwijnen, A.: Measuring the elastic modulus of snow, Geophys. Res. Lett., 44, 11088–11096, https://doi.org/10.1002/2017GL075110, 2017.
Gibson, L. J. and Ashby, M. F.: Cellular Solids: structure and properties, 2nd ed., Cambridge University Press, Cambridge, UK, 510 pp., 1997.
Hagenmuller, P., Theile, T. C., and Schneebeli, M.: Numerical simulation of microstructural damage and tensile strength of snow, Geophys. Res. Lett., 41, 86–89, https://doi.org/10.1002/2013gl058078, 2014.
Hagenmuller, P., Chambon, G., and Naaim, M.: Microstructurebased modeling of snow mechanics: a discrete element approach, The Cryosphere, 9, 1969–1982, https://doi.org/10.5194/tc919692015, 2015.
Johnson, J. B. and Hopkins, M. A.: Identifying microstructural deformation mechanisms in snow using discreteelement modeling, J. Glaciol., 51, 432–442, https://doi.org/10.3189/172756505781829188, 2005.
Kadau, D. and Herrmann, H. J.: Density profiles of loose and collapsed cohesive granular structures generated by ballistic deposition, Phys. Rev. E, 83, 031301, https://doi.org/10.1103/PhysRevE.83.031301, 2011.
KenIchi, K.: Distribution of directional data and fabric tensors, Int. J. Eng. Sci., 22, 149–164, https://doi.org/10.1016/00207225(84)900909, 1984.
Löwe, H., Egli, L., Bartlett, S., Guala, M., and Manes, C.: On the evolution of the snow surface during snowfall, Geophys. Res. Lett., 34, L21507, https://doi.org/10.1029/2007GL031637, 2007.
McClung, D. M.: Shear fracture precipitated by strain softening as a mechanism of dry slab avalanche release, J. Geophys. Res., 84, 3519–3526, https://doi.org/10.1029/JB084iB07p03519, 1979.
Mede, T., Chambon, G., Hagenmuller, P., and Nicot, F.: Snow failure modes under mixed loading, Geophys. Res. Lett., 45, 13351–13358, https://doi.org/10.1029/2018GL080637, 2018.
Mulak, D. and Gaume, J.: Numerical investigation of the mixedmode failure of snow, Comput. Part. Mechan., 6, 439–447, https://doi.org/10.1007/s40571019002245, 2019.
Potyondy, D. O. and Cundall, P. A.: A bondedparticle model for rock, Int. J. Rock Mech. Min., 41, 1329–1364, https://doi.org/10.1016/j.ijrmms.2004.09.011, 2004.
Reiweger, I., Gaume, J., and Schweizer, J.: A new mixedmode failure criterion for weak snowpack layers, Geophys. Res. Lett., 42, 1427–1432, https://doi.org/10.1002/2014GL062780, 2015.
Reuter, B. and Schweizer, J.: Describing snow instability by failure initiation, crack propagation, and slab tensile support, Geophys. Res. Lett., 45, 7019–7027, https://doi.org/10.1029/2018GL078069, 2018.
Reuter, B., Schweizer, J., and van Herwijnen, A.: A processbased approach to estimate point snow instability, The Cryosphere, 9, 837–847, https://doi.org/10.5194/tc98372015, 2015.
Schweizer, J., Jamieson, J. B., and Schneebeli, M.: Snow avalanche formation, Rev. Geophys., 41, 1016, https://doi.org/10.1029/2002RG000123, 2003.
Shertzer, R. H.: Fabric tensors and effective properties of granular materials with application to snow, Ph.D., Montana State University, Bozeman MT, USA, 258 pp., 2011.
Shertzer, R. H., Adams, E. E., and Schneebeli, M.: Anisotropic thermal conductivity model for dry snow, Cold Reg. Sci. Technol., 67, 122–128, https://doi.org/10.1016/j.coldregions.2011.09.005, 2011.
Sigrist, C.: Measurement of fracture mechanical properties of snow and application to dry snow slab avalanche release, Department of Mechanical and Process Engineering, ETH Zurich, Zurich, Switzerland, 139 pp., 2006.
Sigrist, C. and Schweizer, J.: Critical energy release rates of weak snowpack layers determined in field experiments, Geophys. Res. Lett., 34, L03502, https://doi.org/10.1029/2006GL028576, 2007.
Srivastava, P. K., Chandel, C., Mahajan, P., and Pankaj, P.: Prediction of anisotropic elastic properties of snow from its microstructure, Cold Reg. Sci. Technol., 125, 85–100, https://doi.org/10.1016/j.coldregions.2016.02.002, 2016.
van Herwijnen, A. and Jamieson, B.: Highspeed photography of fractures in weak snowpack layers, Cold Reg. Sci. Technol., 43, 71–82, https://doi.org/10.1016/j.coldregions.2005.05.005, 2005.
van Herwijnen, A., Schweizer, J., and Heierli, J.: Measurement of the deformation field associated with fracture propagation in weak snowpack layers, J. Geophys. Res., 115, F03042, https://doi.org/10.1029/2009JF001515, 2010.
van Herwijnen, A., Gaume, J., Bair, E. H., Reuter, B., Birkeland, K. W., and Schweizer, J.: Estimating the effective elastic modulus and specific fracture energy of snowpack layers from field experiments, J. Glaciol., 62, 997–1007, https://doi.org/10.1017/jog.2016.90, 2016.