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

Dry-snow 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 three-dimensional discrete element method (3-D 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 3-D-DEM snow model allowed reproduction of the macroscopic behavior observed during compression and mixed-mode loading of dry-snow 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 load-controlled failure experiments were performed on small samples and compared to results from load-controlled 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) -
Supplement
(485 KB) - BibTeX
- EndNote

Dry-snow 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 so-called 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 pre-defined weak layer until the critical crack length is reached and self-propagation 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 3-D 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 2-D simulations of the PST in agreement with field experiments. However, the oversimplified shape (triangular structure) and the 2-D 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, microstructure-based DEM models adequately reproduce the mechanical behavior (Mede et al., 2018). However, the computational costs of these complex 3-D 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 3-D-DEM 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 high-performance computer (HPC). To relate DEM parameters to macroscopic snow behavior, we will validate the model by simulating basic load cases. Finally, we numerically simulate mixed-mode 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 granular-like 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 parallel-bond contact model (PBM) introduced by
Potyondy and Cundall (2004). The PBM provides the mechanical
behavior of a finite-sized piece of cement-like 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 parallel-bond 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 post-failure phase.

### 2.1.3 System generation

The simulated three-dimensional 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 trade-off 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 snow-like
layers (weak layer or slab layer), unconfined load-controlled 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 Load-controlled test

Load-controlled 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), mixed-mode 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 2nd-order contact tensor coefficients are
defined in Ken-Ichi (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 load-controlled failure experiments on artificially created,
layered snow samples, consisting of a weak layer of depth hoar crystals
between harder layers of fine-grained 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 signal-to-noise
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 mixed-mode 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
load-controlled 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 large-deformation, unconfined load-controlled 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 quasi-brittle 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 load-controlled 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 mixed-mode 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 2nd-order 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 3-D 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 load-controlled compression test, was very similar to that reported by Mede et al. (2018), who simulated snow behavior with microstructure-based snow samples. More generally, Gibson and Ashby (1997) also described these four distinct phases for elastic–brittle foam samples.

The unconfined load-controlled tests under mixed-mode 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 post-failure behavior. Still, as shown above, the post-failure 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 post-failure 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 dry-snow 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 dry-snow slab avalanche.

We developed a mesoscale (between snow grain and slope scale) simulation tool based on 3-D 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 dry-snow 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 mixed-mode 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 mixed-mode loading simulations under different loading angles, was in good agreement with previous experimental and numerical results. The closed-form 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 mixed-mode 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 dry-snow 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/tc-14-39-2020-supplement.

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 co-authors.

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/tc-8-1407-2014, 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/tc-9-1915-2015, 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/tc-11-217-2017, 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/s41467-018-05181-w, 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.: Microstructure-based modeling of snow mechanics: a discrete element approach, The Cryosphere, 9, 1969–1982, https://doi.org/10.5194/tc-9-1969-2015, 2015.

Johnson, J. B. and Hopkins, M. A.: Identifying microstructural deformation mechanisms in snow using discrete-element 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.

Ken-Ichi, K.: Distribution of directional data and fabric tensors, Int. J. Eng. Sci., 22, 149–164, https://doi.org/10.1016/0020-7225(84)90090-9, 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 mixed-mode failure of snow, Comput. Part. Mechan., 6, 439–447, https://doi.org/10.1007/s40571-019-00224-5, 2019.

Potyondy, D. O. and Cundall, P. A.: A bonded-particle 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 mixed-mode 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 process-based approach to estimate point snow instability, The Cryosphere, 9, 837–847, https://doi.org/10.5194/tc-9-837-2015, 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.: High-speed 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.