Microstructure-based modeling of snow mechanics: a discrete element approach

Abstract

Rapid and large deformations of snow are mainly controlled by grain rearrangements, which occur through the failure of cohesive bonds and the creation of new contacts. We exploit a granular description of snow to develop a discrete element model based on the full 3-D microstructure captured by microtomography. The model assumes that snow is composed of rigid grains interacting through localized contacts accounting for cohesion and friction. The geometry of the grains and of the intergranular bonding system are explicitly defined from microtomographic data using geometrical criteria based on curvature and contiguity. Single grains are represented as rigid clumps of spheres. The model is applied to different snow samples subjected to confined compression tests. A detailed sensitivity analysis shows that artifacts introduced by the modeling approach and the influence of numerical parameters are limited compared to variations due to the geometry of the microstructure. The model shows that the compression behavior of snow is mainly controlled by the density of the samples, but that deviations from a pure density parameterization are not insignificant during the first phase of deformation. In particular, the model correctly predicts that, for a given density, faceted crystals are less resistant to compression than rounded grains or decomposed snow. For larger compression strains, no clear differences between snow types are observed.

Introduction

Knowledge of the mechanical properties of snow is required for many applications such as avalanche risk forecasting, optimizing over-snow vehicle traffic, quantifying loads on structures, etc. . Snow is a material composed of air, ice and sometimes liquid water and impurities. Snowflakes can take numerous different forms controlled by the atmospheric conditions e.g.,. Once deposited on the ground, snow continues to evolve. Indeed, its high porosity and its temperature close to its melting point promote rapid changes in the microstructure through water vapor transport or melting/refreezing processes e.g.,. As a consequence, the snow material is characterized by a wide range of densities and microstructural patterns classified into different snow types . In addition, the ice matrix of the snow microstructure exhibits different types of mechanical behavior such as elasticity, visco-plasticity and brittle failure, depending on the applied load, strain rate and temperature . This variety of snow types and deformation mechanisms is a major scientific obstacle to improving knowledge on the mechanical properties of snow .

Numerous studies have attempted to measure the macroscopic mechanical behavior of different snow types under various loading conditions (strain rate and type of experiment) and temperatures. , and reviewed the main results of these experimental studies. Global trends of increasing strength and stiffness with density have been reported. However, plots of mechanical properties versus density generally show a wide scatter . Density alone therefore appears to be insufficient to fully characterize snow's mechanical behavior , although it continues to be used as the main microstructure descriptor in the majority of existing studies. It appears necessary to better account for the microstructure in order to derive more relevant indicators of the mechanical properties of snow .

Models incorporating state variables describing microstructural features such as grain and bond size were developed based on 2-D stereological observations e.g.,. However, these 2-D-based models inevitably introduce a biased representation of the snow microstructure e.g.,. With micro-computed tomography ($\mathit{\mu }$CT), it is now possible to obtain a 3-D image of the microstructure at resolutions down to a few microns . The images obtained can be used as direct input for mechanical simulations . These $\mathit{\mu }$CT-based mechanical models nevertheless tend to be computationally extremely expensive, especially when the microstructure is meshed with finite elements (FE) of similar sizes and no degree of freedom is blocked. Moreover, such FE simulations are inherently limited to relatively small strains and simple material models for ice.

Rapid and large deformations of snow are mainly controlled by grain rearrangements, which are enabled by bond failure and creation. This type of deformation is involved in the release of slab avalanches , in the characterization of a snow profile with an indenter , or when a vehicle wheel comes in contact with snow on the ground . To model this specific deformation regime, we propose to describe snow as a granular material, i.e., to simplify the microstructure into a set of rigid grains interacting only through their contacts. This assumption of snow as a granular material can be considered as reasonable for strain rates greater than ${10}^{-\mathrm{4}}$ s${}^{-\mathrm{1}}$ and for natural seasonal snow with a relatively high porosity ($>0.6$) . It drastically reduces the number of degrees of freedom in the microstructure and makes it possible to simulate snow deformation with the discrete element method (DEM), which is better suited than FE to model large discontinuous deformations governed by microscale processes such as bond breakage and compaction of broken fragments. To our knowledge, the first attempt in this direction was conducted by who investigated snow mechanics with DEM using a simplified representation of the microstructure based on random cylinders with hemispherical ends. Recently, used DEM to study the interaction between snow and a tire tread, snow microstructure being represented as a random assembly of spheres deposited by gravity.

In the present paper, we propose to model dry snow deformation at high loading rates using the DEM approach based on the real 3-D microstructure of snow directly captured by $\mathit{\mu }$CT. In the first part, the different microstructure images used are presented. The DEM model, which requires a method to translate the microstructural information into discrete elements, and the definition of a contact law, is then described. The model capabilities are evaluated on one particular type of loading, namely, uniaxial confined compression, obtained through an imposed uniaxial strain and lateral periodic boundary conditions. This loading involves the creation of new contacts and particle rearrangements, which the model aims to capture. The third section begins with a detailed analysis of the model's sensitivity to its parameters. Lastly, the model is applied to the different snow samples and the simulated compressive behavior is analyzed.

Data set and modeling procedure

3-D image data set

Seven different $\mathit{\mu }$CT snow images were considered, spanning different snow types (decomposing and fragmented precipitation particles (DF), rounded grains (RG), faceted crystals (FC), depth hoar (DH); ). These images were obtained from controlled cold-room experiments (samples s-DFRG, s-RG0, s-RG1 measured by during an isothermal metamorphism experiment; sample s-FC measured by during a temperature–gradient experiment) or field sampling (sample s-DF measured by and samples s-FCDF and s-FCDH measured by ). All samples were scanned with X-ray absorption tomography after impregnation of the material with 1-chloronaphtalene . The $\mathit{\mu }$CT images have a voxel resolution ranging between 5 and 10 $\mathrm{\mu }\mathrm{m}$. The grayscale images were binary segmented with the technique described by with an effective resolution of 2 voxels. The characteristics of the different images are summarized in Table . Typical snow microstructural patterns used in this study are shown in Fig. .

Description of the microtomographic images used in this study. All images are cubic. Density and equivalent spherical radius were computed directly from the binary images. The equivalent spherical radius $r_{\mathrm{eq}}$ is defined as $\mathrm{3}V/S$, with $V$ the volume of ice and $S$ the snow surface area. It represents a typical length scale of the ice matrix.

Sample name	Sample code	Type	Voxel size ($\mathrm{\mu }$m)	Side length (mm)	Density (kg m${}^{-\mathrm{3}}$)	Equivalent spherical radius ($\mathrm{\mu }$m)
s-DF	Ip04	DF	8.59	5.84	145	128
s-DFRG	I08	DF/RG	9.82	4.91	145	111
s-FCDF	4	FC/DF	9.65	5.79	163	191
s-FC	2A	FC	8.37	5.00	319	183
s-FCDH	9	FC/DH	9.65	5.79	183	179
s-RG0	I15	RG	9.82	4.91	182	146
s-RG1	I23	RG	9.82	4.91	242	191

Example of different snow microstructural patterns used in this study: (a) sample s-DF, (b) sample s-FCDH and (c) sample s-RG0. The scale is the same for the three images.

[Figure omitted. See PDF]

Snow microstructure representation

In practice, the output of tomography is a binary image composed of air and a continuous ice matrix, without clearly identifiable grains. The first step thus consists in detecting grains in these binary images. Since the grains are not well separated but rather sintered together, their definition is generally not objective. In the present study, the algorithm developed by was used and is briefly described in the following. In the second step, we explain how the detected grains, which can have complex shapes, are represented in terms of discrete elements by clumping spheres together.

Grain segmentation

developed an algorithm to segment the snow microstructure into grains that are mechanically relevant. The algorithm detects the zones of reduced ice thickness, corresponding to the bonds, based on geometrical criteria on curvature and contiguity. First, potential bond candidates are identified using the value of the minimal principal curvature ${\mathit{\kappa }}_{min}$ in the ice matrix. The bond candidates are located on concavities (${\mathit{\kappa }}_{min}<\mathrm{0}$). A parameter $\stackrel{\mathrm{̃}}{{\mathit{\kappa }}_{\mathrm{t}}}$ sets the sensitivity of this pre-segmentation to curvature: a low value of $\stackrel{\mathrm{̃}}{{\mathit{\kappa }}_{\mathrm{t}}}$ yields numerous bond candidates, whereas a higher value of $\stackrel{\mathrm{̃}}{{\mathit{\kappa }}_{\mathrm{t}}}$ yields fewer bond candidates, which are then located only on the most pronounced concavities. Second, the bond candidates are refined: their area is minimized and highly contiguous grains, i.e., grains that share a large proportion of their total surface, are merged. A contiguity parameter $c_{\mathrm{t}}$ defines the threshold for the merging: a grain is merged to its neighbor if $S_{\mathrm{int}}<c_{\mathrm{t}}\times S_{\mathrm{tot}}$, where $S_{\mathrm{int}}$ is the contact area between the two grains and $S_{\mathrm{tot}}$ is the total grain surface area. showed that this grain segmentation algorithm is able to detect the bonds that support the highest stresses locally. With the assumption of homogeneous ice strength, the algorithm detects the potential failure paths in the microstructure and is thus relevant for the DEM approach.

Clumps of spheres

The segmented grains are represented by sets of voxels. This description must then be translated into discrete elements that can be handled by a DEM model. In principle, clumping of basic geometrical units such as spheres e.g., or considering polyhedral elements e.g., can be relevant solutions to capture the morphology of such irregular 3-D particles. We chose the former approach, namely, the sphere clumping approach, because it is much simpler in terms of numerical implementation. In particular, the clumping approach is free of the singularity issues involved in the definition of contacts between faces, edges and points, which are difficult to handle and computationally demanding in terms of calculation .

In the present study, a very simple technique based on non-overlapping spheres was used to describe (or mesh) the snow grains into discrete elements. Following grain segmentation, the resolution of the 3-D image was reduced by an adjustable factor. Then all voxels at the interface between two grains or between air and ice were replaced by a sphere with a diameter equal to the side length of the voxel with the same center position (Fig. ). The sphere diameter thus depends on the reduction of the initial image resolution. Each grain decomposed into numerous spheres behaves as a rigid body. The contact between grains is computed via the detection of contacts between the spheres. Considering spheres only at the surface of the grains diminishes the computational cost. Owing to the use of equally sized spheres, relatively large time steps can be taken for the simulations.

Granular description of the ice matrix. The binary image (a) is segmented into grains (b) separated by bonds (in black). The grains detected are themselves represented as clumps of spheres (c).

[Figure omitted. See PDF]

Discrete element model

Discrete element modeling was performed with the Yade DEM code (https://yade-dem.org/doc/ ) using a soft contact approach . Each grain is represented by a single clump, itself composed of a rigid aggregate of multiple spheres, called clump members. The motion equations are solved on the clump. The positions of its members are updated according to the clump's overall motion. The forces and moments acting on the clump are the sum of the forces acting on its members. In the model, only contact forces are used and no volumetric forces are considered (no gravity). The interactions between the members of two different clumps are elastic, frictional and cohesive, cohesion being active only for the contacts that exist in the initial undeformed microstructure. Sintering during deformation is therefore not taken into account.

Contact law

In this section, the contact law between two spheres of radius $r$ is detailed, and the resulting contact law between two grains (or clumps) is briefly described. The material, ice, is described by a Young's modulus $E$, a Poisson ratio $\mathit{\nu }$, a normal cohesion ${\mathit{\sigma }}_{\mathrm{n}}$, a shear cohesion ${\mathit{\sigma }}_{\mathrm{t}}$, and a friction angle $\mathit{\phi }$. The force ${\mathbf{f}}_{\mathrm{n}}$ acting on a sphere $S_{\mathrm{1}}$ in contact with a sphere $S_{\mathrm{2}}$ along the normal direction $\mathbf{n}$ pointing from $S_{\mathrm{1}}$ to $S_{\mathrm{2}}$ depends on $E$, ${\mathit{\sigma }}_{\mathrm{n}}$ and the overlap $\mathit{\delta }$ between the spheres as follows: $$\begin{array}{ccc}{\displaystyle }{\mathbf{f}}_{\mathrm{n}}& {\displaystyle }=& {\displaystyle }-\left(Er\mathit{\delta }\right)\mathbf{n}\text{ if }-Er\mathit{\delta }<\mathrm{4}{\mathit{\sigma }}_{\mathrm{n}}{r}^{\mathrm{2}}\\ & {\displaystyle }{\mathbf{f}}_{\mathrm{n}}& {\displaystyle }=& {\displaystyle }\mathbf{0}\text{ and cohesion is broken}\text{ if }-Er\mathit{\delta }\ge \mathrm{4}{\mathit{\sigma }}_{\mathrm{n}}{r}^{\mathrm{2}},\end{array}$$ where $r$ denotes the sphere radius and a positive (resp. negative) sign of $\mathit{\delta }$ corresponds to compression (resp. tension). Equation () simply describes a linear spring with a stiffness $K_{\mathrm{n}}=Er$ and a maximum force $\mathrm{4}{\mathit{\sigma }}_{\mathrm{n}}{r}^{\mathrm{2}}$ in tension before breakage (Fig. a). The force ${\mathbf{f}}_{\mathrm{t}}$ acting on $S_{\mathrm{1}}$ in the tangent plane depends on $E$, $\mathit{\nu }$, ${\mathit{\sigma }}_{\mathrm{t}}$, $\mathit{\phi }$ and the accumulated shear displacement ${\mathbf{u}}_{\mathrm{t}}$ (Fig. b). Let us define the elastic shear force as ${\mathbf{f}}_{\mathrm{t}}^{\mathrm{e}}=-\left(\mathit{\nu }Er\right){\mathbf{u}}_{\mathrm{t}}$ and the frictional shear force as ${\mathbf{f}}_{\mathrm{t}}^{\mathrm{s}}=({\mathbf{f}}_{\mathrm{t}}^{\mathrm{e}}/|{\mathbf{f}}_{\mathrm{t}}^{\mathrm{e}}\left|\right)\left|{\mathbf{f}}_{\mathrm{n}}\right|\tan \left(\mathit{\phi }\right)$. The force ${\mathbf{f}}_{\mathrm{t}}$ is then $$\begin{array}{cc}& {\mathbf{f}}_{\mathrm{t}}={\mathbf{f}}_{\mathrm{t}}^{\mathrm{e}}\text{ if }\left|{\mathbf{f}}_{\mathrm{t}}^{\mathrm{e}}\right|<\mathrm{4}{\mathit{\sigma }}_{\mathrm{t}}{r}^{\mathrm{2}}+\left|{\mathbf{f}}_{\mathrm{n}}\right|\tan \left(\mathit{\phi }\right)\\ & & {\mathbf{f}}_{\mathrm{t}}={\mathbf{f}}_{\mathrm{t}}^{\mathrm{s}}\text{ and cohesion is broken}\text{ if }\left|{\mathbf{f}}_{\mathrm{t}}^{\mathrm{e}}\right|\ge \mathrm{4}{\mathit{\sigma }}_{\mathrm{t}}{r}^{\mathrm{2}}+\left|{\mathbf{f}}_{\mathrm{n}}\right|\tan \left(\mathit{\phi }\right).\end{array}$$ The tangential force ${\mathbf{f}}_{\mathrm{t}}$ also adds a torque ${\mathbf{t}}_{\mathrm{t}}=d(-\mathbf{n})\times {\mathbf{f}}_{\mathrm{t}}$ on sphere $S_{\mathrm{1}}$, with $d=r-\mathit{\delta }/\mathrm{2}$ the distance from the center of $S_{\mathrm{1}}$ to the contact point. Note that if cohesion is broken either in tension or in shear, then both ${\mathit{\sigma }}_{\mathrm{t}}$ and ${\mathit{\sigma }}_{\mathrm{n}}$ are set to zero on this contact. The new contacts have no cohesion.

In this study, cohesion strengths in tension (${\mathit{\sigma }}_{\mathrm{n}}$) and in shear (${\mathit{\sigma }}_{\mathrm{t}}$) are assumed to be equal and are simply denoted ${\mathit{\sigma }}_{\mathrm{ice}}$ (${\mathit{\sigma }}_{\mathrm{n}}={\mathit{\sigma }}_{\mathrm{t}}={\mathit{\sigma }}_{\mathrm{ice}}$). Furthermore, we use a constant Poisson ratio of 0.3. Hence the considered contact law is fully parameterized by the contact stiffness $E$ (elastic part), the strength ${\mathit{\sigma }}_{\mathrm{ice}}$ (cohesive part) and the friction coefficient $\tan \mathit{\phi }$ linking the shear force resisting to sliding to the applied normal force (frictional part).

Considering now intergranular contacts, they are generally composed of several sphere–sphere contacts as described above (Fig. c and d). For normal loading and simple shearing, the overall cohesion force of a granular contact of area $A$ is the sum of the cohesion force $\mathrm{4}{\mathit{\sigma }}_{\mathrm{ice}}{r}^{\mathrm{2}}$ of the $N=A/\left(\mathrm{4}{r}^{\mathrm{2}}\right)$ equally loaded sphere–sphere contacts (Fig. c). This overall cohesion is thus equal to ${\mathit{\sigma }}_{\mathrm{ice}}A$. For friction, the bumpy surface at the contacts between grains yields an effective granular friction coefficient varying with displacement (Fig. d). However, for the parameters considered in this study, the mean effective intergranular friction coefficient remains close to the microscopic friction coefficient defined at the sphere–sphere interactions ($\tan \left(\mathit{\phi }\right)$).

Behavior of a sphere–sphere contact (a, b) and a grain–grain contact (c, d) in two simple cases: normal loading (a, c) and simple shearing (b, d) under constant normal force (${\mathbf{f}}_{\mathrm{n}}$, ${\mathbf{F}}_{\mathrm{n}}$). The grain–grain contact is here composed of seven sphere–sphere contacts. The small spheres have a radius $r$. The intergranular contact area is denoted $A$. The number of sphere–sphere interactions in the grain–grain contact is denoted $N=A/\left(\mathrm{4}{r}^{\mathrm{2}}\right)$.

[Figure omitted. See PDF]

Motion equation integration

A sufficiently short time step d$t$ is necessary for the stability of the algorithm and the physical relevance of the simulation. Usually the critical time step d$t_{\mathrm{c}}$ is estimated from the propagation speed of elastic waves in the sample, on the basis of a single degree-of-freedom mass $m$ connected to ground by a linear spring of stiffness $k$: d$t_{\mathrm{c}}\approx \sqrt{m/k}$ . Here, $m$ corresponds to the smallest grain mass, which is at least composed of one sphere, and $k$ is the spring coefficient of a sphere–sphere interaction. Hence, the critical time step can be estimated as $$\mathrm{d}{t}_{\mathrm{c}}\approx \sqrt{{\mathit{\rho }}_{\mathrm{ice}}{r}^{\mathrm{3}}/\left(Er\right)}=r\sqrt{{\mathit{\rho }}_{\mathrm{ice}}/E}.$$ Since the main focus of this study is on deformation regimes controlled by grain rearrangements, the elasticity of the contacts is expected to have little influence on the macroscopic response of the sample in the domain of interest e.g.,. Hence, the value of $E$ was decreased compared to the real stiffness of ice, and the value of ${\mathit{\rho }}_{\mathrm{ice}}$ was increased, in order to increase d$t_{\mathrm{c}}$ and speed up the simulations. The contact stiffness $E$ must nevertheless remain large enough to ensure that overlaps remain small compared to the sphere radius (assumption of rigid grains). The maximum stress at the contacts is expected to be on the order of the cohesion (${\mathit{\sigma }}_{\mathrm{ice}}$), which yields an overlap $\mathit{\delta }\approx r\cdot (E/{\mathit{\sigma }}_{\mathrm{ice}})$. Accordingly, a value $E=10\cdot {\mathit{\sigma }}_{\mathrm{ice}}$, corresponding to $\mathit{\delta }/r<0.1$, was chosen to expedite the simulations. As will be discussed later (Sect. ), this value of $E$ is slightly too low to completely fulfill the rigid grain assumption, but the effect of $E$ was considered sufficiently small, at least in a first step. Similarly, the density of ice ${\mathit{\rho }}_{\mathrm{ice}}$ was taken as ${10}^{\mathrm{4}}$ kg m${}^{-\mathrm{3}}$ (instead of 917 kg m${}^{-\mathrm{3}}$ for real ice). We checked that this artificial adjustment of the grain mass does not affect the overall simulated mechanical behavior, since the inertial effects are negligible under the loading conditions considered (Sect. ). These choices allowed us to use a time step d$t=\mathrm{5}\times {10}^{-\mathrm{7}}s<\mathrm{d}{t}_{\mathrm{c}}$.

To dissipate kinetic energy and avoid numerical instabilities, a numerical damping (coefficient ${\mathit{\lambda }}_a)$ was also used . The basic idea of this damping is to decrease forces that increase particle velocities, and vice versa. A value ${\mathit{\lambda }}_a=0.2$ was used in this study. In practice, this parameter appears to have little effect on the simulated mechanical behavior as long as the simulations are quasi-static.

The simulations presented in this paper were run on a desktop computer with a single processor (2.7 GHz) and 16 Gb RAM. The typical computing time of the simulation is on the order of 1 day. More precisely, the test shown in Sect.  runs on this computer in 15 h with, in average, 8 time-step iterations per second. It requires 2.5 Gb RAM. The sample s-RG0 of rounded grains (${\mathrm{4}}^{\mathrm{3}}$ mm${}^{\mathrm{3}}$) is composed of 800 grains, and described with about ${10}^{\mathrm{5}}$ spheres ($1.5\times {10}^{\mathrm{5}}$, if the spheres located inside the grains are not removed). On this sample, the mean number of sphere–sphere contacts per cohesive intergranular bond is 14.

Boundary conditions

The boundary conditions were chosen to reproduce strain-controlled vertical confined compression of snow. Hence, the samples were placed between two horizontal plates discretized using spheres of the same size as the spheres used to describe the grain shape (Fig. ). The bottom plate was fixed and a constant vertical velocity $V_{\mathrm{plate}}$ was applied to the top plate. The contact law between the grains and the plates is the same as the intergranular contact law. In the horizontal directions, periodic boundary conditions were applied.

The velocity of the upper plate $V_{\mathrm{plate}}$ prescribes the compressive strain rate imposed on the sample. A value of $V_{\mathrm{plate}}={10}^{-\mathrm{2}}$ m s${}^{-\mathrm{1}}$ was chosen, checking that kinetic energy always remained negligible compared to the total energy in the sample (see Sect. ). Therefore, the simulated mechanical behavior can be regarded as quasi-static, and independent of $V_{\mathrm{plate}}$. In addition, let us recall that the granular-like behavior of snow is valid only for relatively high strain rates ($\dot{\mathit{ϵ}}>{10}^{-\mathrm{4}}$ s${}^{-\mathrm{1}}$, ), a condition that is largely fulfilled with the value of $V_{\mathrm{plate}}$ considered.

Boundary conditions used in the present study. The grains (medium gray) are compressed between two plates composed of spheres (black): the bottom plate is fixed and the top plate is moving down with a constant velocity $V_{\mathrm{plate}}$. Periodic space was used in the horizontal directions. Elements (hatched) crossing the periodic cell limits (dotted lines) appear on the opposite side of the samples. Note that this figure represents a 2-D simulation for clarity purposes, but the simulations performed are 3-D.

[Figure omitted. See PDF]

Results

First, a typical simulation performed with the sample of rounded grains (s-RG1) is described to give an overview of the deformation mechanisms occurring during compression loading. A sensitivity analysis of the model to the DEM parameters and microstructure representation is then presented. Lastly, the model is applied to all snow images described in Table .

Deformation phases

The model parameters used for this first simulation are summarized in Table . The sample undergoes uniaxial compression with an imposed global vertical strain ${\mathcal{E}}_{zz}$ defined as the ratio of the displacement of the upper plate to the initial height $L_{\mathrm{0}}$ of the sample. In response, the sample produces a vertical reaction force $F_z$, which defines an overall vertical stress ${\mathrm{\Sigma }}_{zz}=F_z/L_{\mathrm{0}}^{\mathrm{2}}$, where $L_{\mathrm{0}}^{\mathrm{2}}$ is the horizontal cross sectional area of the sample. The relation between the imposed macroscopic strain ${\mathcal{E}}_{zz}$ and the resulting macroscopic stress ${\mathrm{\Sigma }}_{zz}$ is shown in Fig. a. When the sample is compressed, its density $\mathit{\rho }$ increases as $\mathit{\rho }={\mathit{\rho }}_{\mathrm{0}}/(\mathrm{1}-{\mathcal{E}}_{zz})$, where ${\mathit{\rho }}_{\mathrm{0}}$ is the initial density. The work of the imposed deformation adds mechanical energy to the sample, which is stored as kinetic and elastic energy, and dissipated through plasticity (shear sliding), bond breaking and non-viscous damping (Fig. b). During deformation, cohesive bonds are broken and new non-cohesive contacts are created (Fig. c). The distributions of overlap $\mathit{\delta }$ between spheres are shown for different strain values on Fig. d. The mechanisms observed can be classified into three phases corresponding to different macroscopic strain intervals.

Model parameters used for the simulations presented in this paper.

Name	Symbol	Value
Material properties
Young's modulus	$E$	10 MPa
Poisson ratio	$\mathit{\nu }$	0.3
Friction angle	$\tan \left(\mathit{\phi }\right)$	0.2
Normal cohesion strength	${\mathit{\sigma }}_{\mathrm{n}}={\mathit{\sigma }}_{\mathrm{ice}}$	1 MPa
Shear cohesion strength	${\mathit{\sigma }}_{\mathrm{t}}={\mathit{\sigma }}_{\mathrm{ice}}$	1 MPa
Density of the grains	${\mathit{\rho }}_{\mathrm{ice}}$	10${}^{\mathrm{4}}$ kg m${}^{-\mathrm{3}}$
Boundary conditions
Side length of the cubic snow sample	$L_{\mathrm{0}}$	$\mathrm{4}$ mm
Compression velocity	$V_{\mathrm{plate}}$	$10$ mm s${}^{-\mathrm{1}}$
Integration of motion equation
Time step	d$t$	$\mathrm{5}\times {10}^{-\mathrm{7}}$ s
Non-viscous damping coefficient	${\mathit{\lambda }}_a$	0.2
Representation of the microstructure
Segmentation parameter on curvature	$\stackrel{\mathrm{̃}}{{\mathit{\kappa }}_{\mathrm{t}}}$	1.0
Segmentation parameter on contiguity	$c_{\mathrm{t}}$	0.25
Radius of the spherical members	$r$	25 $\mathrm{\mu }\mathrm{m}$

Elastic phase for ${\mathcal{E}}_{zz}\mathit{\lesssim }0.02$

After the initiation of the compression, the snow structure first undergoes elastic deformation. The work of the loading is completely stored as elastic energy at the bonds. No cohesive bonds are broken, nor are any new bonds created. If the upper plate is moved back to its initial position, the stress decreases to zero and the microstructure fully recovers its initial state with no residual deformation or stress. As a consequence, a linear relationship between stress ${\mathrm{\Sigma }}_{zz}$ and strain ${\mathcal{E}}_{zz}$ is observed in this elastic phase. The slope of the stress–strain relation quantifies the effective elastic modulus $E_{\mathrm{s}}$ of the granular assembly. It is noteworthy that a significant number of bonds are in tension (negative overlaps exist, as shown in Fig. d), even if the macroscopic loading is compression.

Brittle/frictional phase for $0.02\mathit{\lesssim }{\mathcal{E}}_{zz}\mathit{\lesssim }0.3$

When the strain increases (${\mathcal{E}}_{zz}\mathit{\gtrsim }0.02$), the first bonds start to break and the snow sample deforms with an approximately constant resistance of about 10 kPa. This apparent plastic behavior is caused by the geometrical rearrangements of grains made possible by bond breakage and grain sliding. If the upper plate is moved back to its initial position, permanent deformations are observed. The breaking of bonds leads to small negative jumps of the stress, which are balanced by friction on the broken bonds, the activation of other cohesive bonds and the creation of new bonds, thus resulting in the almost constant macroscopic stress. The friction between grains yields plastic dissipation. Kinetic energy always remains negligible compared to total energy. The number of cohesive contacts appears to strongly decrease during this phase, but the total number of interactions only decreases slightly. The structure tightens because bond failure permits grains to move into the pore space.

Dense compaction phase for $0.3\mathit{\lesssim }{\mathcal{E}}_{zz}$

In this phase, hardening of the structure is observed. The stress progressively increases, following an exponential-like trend. This phase starts when the grain packing has reached a certain density (here $\approx 450$ kg m${}^{-\mathrm{3}}$) and can no longer accommodate the imposed displacement with little resistance. The stress increase is accompanied by the creation of numerous contacts. Because most contacts are not cohesive in this phase, the distribution of overlaps shows a clear majority of bonds loaded in compression. Kinetic energy remains negligible compared to total energy.

Results of the model applied to sample s-RG1 (rounded grains) with the parameters listed in Table : (a) macroscopic stress–strain curve, (b) evolution of energy components with strain, (c) total number of bonds and number of cohesive bonds as a function of strain and (d) occurrence distribution of the relative grain overlap for three different imposed strains.

[Figure omitted. See PDF]

Sensitivity analysis

A specific sensitivity analysis was conducted for the following model parameters (Table ):

• The parameters related to the microstructure representation, namely, the number of grains used to discretize the microstructure (controlled by the parameters ${\stackrel{\mathrm{̃}}{\mathit{\kappa }}}_{\mathrm{t}}$ and $c_{\mathrm{t}}$ of the grain segmentation algorithm) and the size of the spheres used to describe the grain shape (controlled by $r$).

• The parameters related to the material properties, namely, the values of cohesion ${\mathit{\sigma }}_{\mathrm{ice}}$ and microscopic friction coefficient $\tan \left(\mathit{\phi }\right)$. As discussed above, the other parameters (${\mathit{\rho }}_{\mathrm{ice}}$, $E$) appearing in the contact law are expected to have little influence on the results in the brittle/frictional phase (see Sect. ). However, since the values of these parameters were chosen artificially to expedite the computations, a sensitivity study has also been conducted.

• The side length of the sample $L_{\mathrm{0}}$. The model is expected to characterize the mechanical behavior of the sample regardless of its size. For this purpose, the sample needs to be large enough and at least larger than a representative elementary volume (REV). Investigations of the REV size related to the compressive behavior of snow will be presented.

Microstructure representation

Example of the discrete element representation of two grains but with different effective resolutions (diameters of spheres).

[Figure omitted. See PDF]

Simulated stress–strain curves for different sphere sizes (a) and number of grains (b). In (b), the images with 397, 668, 857 and 1132 grains were obtained with the segmentation parameter $({\mathit{\kappa }}_{\mathrm{t}},c_{\mathrm{t}})$ set to $(1.3,0.08)$, $(1.3,0.16)$, $(1.0,0.25)$ and $(0.7,0.5)$, respectively.

[Figure omitted. See PDF]

Figure a shows the stress–strain curves obtained for different sphere sizes, i.e., different reductions of the resolution of the initial 3-D grain image. The simulated mechanical behavior appears relatively insensitive to the sphere size in the range $[20,40]$ $\mathrm{\mu }$m. The stress–strain curve, however, starts to deviate for the largest spheres tested (60 $\mathrm{\mu }$m), due to an insufficiently fine description of the microstructure. Accordingly, a sphere size in the range $[20,40]$ $\mathrm{\mu }$m appears, here, as a reasonable choice to describe the mechanical behavior.

The total number of grains in the sample is not absolute, but depends on the chosen segmentation parameters ${\mathit{\kappa }}_{\mathrm{t}}$ and $c_{\mathrm{t}}$. The same image of the sample s-RG1 was segmented into more or fewer grains with different sets of parameters ${\mathit{\kappa }}_{\mathrm{t}}$ and $c_{\mathrm{t}}$. As shown in Fig. b, the simulated mechanical behavior is not very sensitive to the number of grains in the range $[668,1132]$, but the stress–strain curve obtained with a very low number of grains (397) differs significantly. This reflects the requirement that the chosen decomposition into grains should be sufficiently fine to preserve the microstructure's essential degrees of freedom.

Friction and cohesion at intergranular contacts

reported a static ice–ice coefficient of friction equal to 0.1. Note, however, that the uncertainty on this quantity remains high, and this reported value of $\tan \left(\mathit{\phi }\right)$ only provides an order of magnitude. Figure a shows the simulated stress–strain curve for different values of the friction coefficient between 0.1 and 0.4. It turns out that the value of $\tan \left(\mathit{\phi }\right)$ has little effect on the computed mechanical behavior for macroscopic strains in the range $[\mathrm{0},0.2]$. In contrast, for larger strains, the stress increases significantly with the value of $\tan \left(\mathit{\phi }\right)$. In the so-called dense compaction phase, the increase of non-cohesive intergranular contacts (see Fig. c) thus enhances the effect of the microscopic friction coefficient on the macroscopic mechanical behavior.

reported tensile strength of ice in the range $[0.7,3.1]$ MPa. The data are widely scattered due to the effects of temperature, strain rate, crystal size, etc. As shown in Fig. b, the modeled macroscopic behavior of snow strongly depends on the microscopic cohesion at the contacts. For a given strain ${\mathcal{E}}_{zz}$, an almost linear relationship between ${\mathrm{\Sigma }}_{zz}$ and ${\mathit{\sigma }}_{\mathrm{ice}}$ is observed. In the case of pure tension without creation of new contacts, this linear relation between the microscopic and macroscopic strengths would be obvious. In the present case of compression, however, which involves spatial grain rearrangements and creation of new contacts, this linear relation was not necessarily expected. In the following, we fixed the microscopic cohesion to 1 MPa. In future studies, effects of temperature on the macroscopic behavior of snow, to the extent that sintering effects remain limited, could also be considered by accounting for the influence of temperature on the microscopic cohesion.

Simulated stress–strain curves for different microscopic friction coefficients (a) and cohesion strengths (b). Note that, in (b) the macroscopic stress ($y$ axis) is scaled by the microscopic strength ${\mathit{\sigma }}_{\mathrm{ice}}$ and that the contact stiffness $E$ varies here as $E=10\times {\mathit{\sigma }}_{\mathrm{ice}}$ in order to keep a constant microscopic strain at failure.

[Figure omitted. See PDF]

Contact stiffness and ice density

With the value $E=10$ MPa used in this study, the typical granular interpenetrations encountered during the brittle/frictional phase (Fig. d) are too large (around a few percents of the grain radius) to completely ensure the assumption of rigid grains. Accordingly, Fig. a shows that the mechanical response of the sample is not completely independent of $E$ in this phase. However, this sensitivity on $E$ remains limited compared to the typical “noise” level on the stress–strain curves, and was considered as acceptable in this study. We note that simulations with $E=100$ MPa take as much as 1 week computational time compared to simulations with $E=10$ MPa, which becomes prohibitive when trying to compare the response of different snow types. As expected, we also observe that the influence of $E$ dramatically increases in the dense compaction phase.

In order to expedite the simulations, the density of ice ${\mathit{\rho }}_{\mathrm{ice}}$ was, here, taken as ${10}^{\mathrm{4}}$ kg m${}^{-\mathrm{3}}$ instead of 917 kg m${}^{-\mathrm{3}}$. Figure b shows the sensitivity of the computed mechanical behavior for different values of ${\mathit{\rho }}_{\mathrm{ice}}$. As shown, this artificial adjustment of ${\mathit{\rho }}_{\mathrm{ice}}$ to ${10}^{\mathrm{4}}$ kg m${}^{-\mathrm{3}}$ does not affect the overall simulated mechanical behavior. This is explained by the very low inertial effects on the overall mechanical behavior.

Computed stress–strain curves for different Young's moduli (a) and ice densities (b).

[Figure omitted. See PDF]

Representative elementary volume

The REV can be defined as the smallest fraction of the sample volume over which the measurement or simulation of a given variable will yield a value representative of the macroscopic behavior. We assessed the representativity of mechanical simulations performed on volumes of about $\mathrm{5}\times \mathrm{5}\times \mathrm{5}$ mm${}^{\mathrm{3}}$, by considering the results obtained with smaller sub-volumes. Two snow samples, namely, s-DF and s-RG1, which span the tested range of snow density (Table ), were considered. As shown in Fig. , the simulated stress–strain curves progressively converge with increasing volume. For the smallest tested volumes, stress fluctuations associated with the failure of single bonds are so large that it is not possible to define a homogeneous mechanical behavior. For volumes larger than $\mathrm{4}\times \mathrm{4}\times \mathrm{4}$ mm${}^{\mathrm{3}}$, however, the variations of the stress–strain curves with the volume size become negligible compared to the variability associated with the choice of the model parameters. The sample sizes used in this study (see Table ) can therefore be regarded as large enough to be representative of the mechanical behavior of snow under compression. Note that these REV sizes related to compression appear much larger than the REV sizes related to density measurements (${1.5}^{\mathrm{3}}$ mm${}^{\mathrm{3}}$, ), but are in good agreement with the REV sizes associated with geometrical parameters characterizing the microstructural bonding system (${\mathrm{3}}^{\mathrm{3}}-{\mathrm{6}}^{\mathrm{3}}$ mm${}^{\mathrm{3}}$, ).

Simulated stress–strain curves for different volume sizes on sample s-RG1 (a) and sample s-DF (b).

[Figure omitted. See PDF]

Behavior of different snow microstructures

The numerical compression experiments were applied to different snow samples spanning different snow types (see Table ). The same model parameters, summarized in Table , were used for all images. The sample sizes are indicated in Table . The only parameter differing between the samples is the sphere size used (or effective resolution of the 3-D image). Samples s-FCDF, s-FC, s-FCDH and s-RG1 are characterized by similar values of the equivalent spherical radius $r_{\mathrm{eq}}$, which roughly characterizes the mean thickness of the ice matrix (Table ). The sensitivity analysis on sample s-RG1 showed that the representation of the ice matrix with spheres with a 25 $\mathrm{\mu }$m radius, in practice corresponding to a reduction of the 3-D image resolution by a factor of 5, is satisfactory. This value $r=25$ $\mathrm{\mu }\mathrm{m}$ was therefore used for samples s-FCDF, s-FC, s-FCDH and s-RG1. Samples s-DF, s-DFRG and s-RG0 present a finer ice matrix characterized by lower values of the equivalent spherical radius (Table ). These samples were therefore decomposed into spheres with a 20 $\mathrm{\mu }\mathrm{m}$ radius, corresponding to a reduction of the 3-D image resolution by a factor of 4 only.

Figure a shows the simulated stress–strain curves computed for all the samples. The three different phases identified on sample s-RG1 can be distinguished in all cases: elastic phase, brittle/frictional phase and dense compaction phase. After the elastic phase, the first bonds start to break for a macroscopic stress ${\mathrm{\Sigma }}_{zz}$ ranging between 1 and 12 kPa, depending on the snow sample. For sample s-FC, the brittle/frictional phase is almost absent and the structure almost immediately exhibits pronounced hardening. In contrast, samples s-DF, s-DFRG, s-FCDF, s-FCDH and s-RG0 present the plastic brittle/frictional phase and can undergo very large deformations (up to 0.4) with an almost constant stress. As discussed below, the level of this stress plateau appears to slightly differ between the samples. Note also that the results obtained for samples s-DF and s-DFRG, which are both mainly composed of decomposed and fragmented snow with a density of 145 kg m${}^{-\mathrm{3}}$, are very similar, which confirms the consistency of the model.

A remarkable result, shown in Fig. b, is that, when plotted as stress–density relations, all the data appear to collapse, at least as a first approximation, onto a single trend. This result shows that the resistance of snow to compression is mainly a function of density, even if the microstructures tested span very different snow types. Some variability, which cannot be explained by density only, is nevertheless observed between the samples, especially during the brittle/frictional phase. Hence, for similar density values below 300 kg m${}^{-\mathrm{3}}$, samples s-FCDF and s-FCDH, which are mainly composed of faceted crystals, present systematically lower resistances than samples s-DF, s-DFRG and s-RG0, composed of decomposed and fragmented snow or rounded grains. Sample s-RG1, which results from isothermal metamorphism of precipitation particles, presents the highest resistance to compression in the density range $[250,300]$ kg m${}^{-\mathrm{3}}$. The model predictions are therefore consistent with the commonly accepted idea that snow types resulting from a temperature–gradient metamorphism are less resistant to mechanical loading than snow types resulting from an isothermal metamorphism e.g.,.

For density values greater than 300 kg m${}^{-\mathrm{3}}$, no significant trend between the different snow types was recognized. The shape and size of the grains in the initial microstructure do not appear to be determinant during this high compaction phase of the samples. In this phase, the material is largely broken into individual grains so that the mechanical response of the assembly is largely independent of the initial microstructure.

Simulated stress–strain curves (a) and stress–density relation (b) for the different snow samples (Table ). The legend is the same for (a) and (b).

[Figure omitted. See PDF]

Discussion

Identification of potential model artifacts

We discuss here the potential artifacts introduced by the model in the simulation of the mechanical behavior of snow under compression.

Elastic phase

Let us recall that the model is composed of rigid grains that cannot deform. Only the contacts between members of different grains can deform to account for intergranular strain. With FE simulations showed that even if the highest strains (or stresses) are located on the bonds, the overall elastic strain is mainly due to the deformation of the grains themselves, since they represent the vast majority of the ice volume. Therefore, even if a realistic value was used for the Young's modulus of ice at the contacts, we would not expect the present DEM model to provide a correct macroscopic Young's modulus. The objective of the model is not to reproduce the elastic properties of snow, but to capture large deformations governed by grain rearrangements. Furthermore, since the Young's modulus at contacts was chosen relatively small compared to the stiffness of ice (see Sect. ), the strain interval for the elastic phase is overestimated. In real snow, irreversible damage starts to occur for smaller strains of about ${10}^{-\mathrm{4}}$ .

Brittle/frictional phase

This phase is controlled by the progressive failure of cohesive bonds and the friction between grains. Failures of intergranular bonds are determined by the force distribution in the sphere–sphere interactions. According to the sensitivity analysis, the model appears to be insensitive to the size of the spheres used to discretize the shape of the grains, and to the number of grains (in reasonable ranges). The friction between grains depends on the roughness of the grain surfaces and on the microscopic friction coefficient. The regular mesh used in the present study tends to create “bumpy” grain surfaces (Fig. ), but, according to the sensitivity analysis, the size of these bumps does not alter the model's results. Moreover, the model exhibits low sensitivity to the friction coefficient in this phase. These results indicate that both the granular description of the snow microstructure and the representation of grains as clumps of spheres are satisfactory to reproduce the mechanical behavior of the material during the brittle/frictional phase.

Dense compaction phase

In this phase, rearrangements of grains can no longer be accommodated within the pore space. The grains form a dense packing. Above a certain macroscopic stress level, we can expect high intra-granular stresses to develop, which could lead to ice grains breaking into smaller parts. The assumption of “unbreakable” grains thus becomes questionable. Moreover, the high intergranular stresses yield large overlaps between spheres compared to their radius. The assumption of rigid grains is no longer valid and the simulated behavior becomes very sensitive to the value chosen for the Young's modulus.

In summary, the approach retained in this study, consisting in modeling snow as a granular material, is certainly reasonable and promising when attempting to reproduce the brittle/frictional deformation phase but is inappropriate to simulate pure elasticity or high compaction phases.

Influence of density and microstructure on the mechanical behavior of snow

We discuss here the first-order role played by density on the simulated mechanical behavior. As explained in the introduction, density is often described as an insufficient mechanical indicator, because of the substantial scatter observed in property–density relations derived from direct experimental measurements e.g.,. In contrast, in this study, density appears to be a very good indicator of the resistance of snow under compression, with only a little scatter attributable to the snow type. In other microstructure-based models, a similar strong dependence on density was also observed for thermal conductivity , tensile strength or the Young's modulus . The main difference between the microstructure-based simulations and direct experimental measurements is the size of the volume tested. Hence, we argue that the spatial variability of density in the large samples used in experiments might constitute an explanation for the scatter observed in property–density relations, especially for properties that are extremely sensitive to density.

Furthermore, the dominant role played by density in the mechanical behavior simulated is probably also relative to the loading conditions considered. In compression, both the mechanical properties and density are expected to depend on the grain bonding system, thus promoting the existence of an apparent relation between stress and density . The model appears capable of reproducing this feature, and we expect the dependence on density to be less pronounced in the case of different loading conditions, such as shear, for which microstructure properties and anisotropy are expected to be greater.

Conclusions

This paper presents a novel approach to modeling the mechanical behavior of snow under large and rapid deformations based on the complete 3-D microstructure of the material. The geometry of the microstructure is directly translated into discrete elements by accounting for the shape of the grains and the initial bonding system of the ice matrix. The grains are rigid and the overall deformations are only due to the geometric rearrangements of grains made possible by bond failure. The sensitivity analysis of the model to its parameters showed that the effects of the microstructure geometry on the simulated mechanical response are not shadowed by numerical artifacts.

In this study, the model was used to reproduce the mechanical behavior of snow under confined compression. The representative volume related to this type of loading was estimated to be ${\mathrm{4}}^{\mathrm{3}}$–${\mathrm{5}}^{\mathrm{3}}$ mm${}^{\mathrm{3}}$. The numerical results showed that, for this type of loading, the mechanical behavior of snow is mostly controlled by density. The stress–density relationships computed for different samples follow a single trend with little variability, even if the tested microstructures span very different snow types. Nevertheless, noticeable second-order effects of the microstructure are observed during the brittle/frictional deformation phase (past the initial elastic regime). In particular, for a given density, snow samples resulting from temperature–gradient metamorphism appear to be less resistant to compression than those resulting from isothermal metamorphism. The model appears capable of accounting for the role of the microstructure on the mechanical properties. During the latter high compaction phase observed for high-density values, the model did not reveal any clear influence of the grain shape and size.

To further investigate the influence of microstructure on snow mechanical properties, this DEM model could straightforwardly be applied to other loading conditions. In particular, investigating shear loading can be useful in the context of slope stability and avalanche release. Another promising prospect would consist in numerically reproducing the penetration of an indenter in the snow microstructure, with the objective of better understanding the signals delivered by micro-penetrometers, which are routinely used for snow characterization . Lastly, directly conducting mechanical experiments in the $\mathit{\mu }$CT would also be useful to evaluate the model and possibly reduce the uncertainties on the microscopic contact parameters. A large scatter exists in the literature on the values of ice cohesion and friction coefficient, and systematic comparisons between experimental results and numerical simulations could help improving the knowledge on these mechanical parameters.

Acknowledgements

Funding by the VOR research network (Tomo_FL project) is acknowledged. We thank T. Mede for his preliminary work on the model and the scientists of 3SR Laboratory and ESRF ID19 beamline, where the 3-D images were obtained. We also thank the CEN staff, especially F. Flin and B. Lesaffre, for the data acquisition. Irstea and CNRM-GAME/CEN are part of Labex OSUG@2020 (Investissements d'Avenir, grant agreement ANR-10-LABX-0056). Irstea is member of Labex TEC21 (Investissements d'Avenir, grant agreement ANR-11-LABX-0030). We thank M. Hopkins and two anonymous reviewers for their constructive feedback on this work. Edited by: M. Schneebeli