Mechanical behaviour of weak snow layers : modelling a porous structure of sintered grains

Weak layers of snow are thin layers of low density and cohesion that consist of a complex network of sintered ice grains. These layers are very fragile and their collapse is considered to be the primary cause of dry slab snow avalanche release. This study reports on an attempt at modelling the mechanical response of these weak snow layers with discrete element (DEM) simulations, using x-ray tomographical images of real snow samples as input data on the microstructure of the material. An original method for representing irregular grain shapes in the DEM is first introduced. The method is based on utilizing the medial axis concept to represent an arbitrary grain shape with an optimal set of overlapping spheres. A thorough study of the effect of the grain approximating technique parameters on the key geometrical features of the grains is then carried out. Finally, the functionality of the model is demonstrated by performing an oedometric test on an image of a real snow


Snow as a granular material
Weak snow layers are thin layers of low cohesion and density, naturally occurring in the snowpack.The understanding of their mechanical behaviour is considered to be of great importance as they are believed to have the decisive effect on a dry slab avalanche release.These layers can be viewed as porous structures of sintered ice grains, that collapse under loading and disintegrate into an cohesionless assembly of grains, withdrawing support to the overlaying slab layer of snow.
Field observations show that weak layers are usually formed by very distinct snow grain types [1].While it is clear that the mechanical character of the weak layer is defined by its microstructure i.e. shape and geometrical arrangement of the snow grains, the precise relation between these properties has not yet been quantified.The extremely fragile nature of these layers makes systematic laboratory experiments notoriously difficult [2,3].
We therefore propose studying the effect of the microstructure on the mechanical behaviour of snow by conducting numerical experiments on a numerical model of a weak layer.The model, naturally has to account for the actual microstructure of the weak layer and has to able to follow the material through different stages -from a porous network of sintered grains to the final stage of a decomposed granular material.
The discrete element method (DEM) is an efficient tool for exploring the mechanics of granular materials.Johnson and Hopkins [4] were the first to demonstrate the e-mail: tijan.mede@irstea.frusefulness of DEM in studying the mechanics of snow.Their simulation was however conducted on a granular deposition of idealized grain shapes and did not account for the actual microstructure of snow.
Recently Hagenmuller [5][6][7] managed to combine DEM approach with x-ray tomography images of snow in order to come up with a method which could account for the microstructure of the material, while having the potential to follow the large deformations and rearrangements within the snow matrix.To this end the tomographical image of snow is divided into individual grains by detecting weak mechanical points i.e. neck regions in the snow matrix.These grains are modelled as the basic elements of the simulation, as separate unbreakable entities, which are sintered into the original matrix of intact snow.When applying mechanical loading to the snow structure, the sintered bonds break and the ice matrix gradually transforms into a cohesionless granular medium.The present study is based on this particular approach to modelling the weak layer of snow, where a novel type of grain geometry representation has been added to improve the computing efficiency of the model.

Approximation of grain geometry
Since the objective of the study is the quantification of the relation between microstructure and the overall mechanical response of the medium, the particular representation of the grains in the DEM simulation is of paramount importance.A grain approximating method is needed that allows on one hand an accurate representation of grain shapes, but on the other hand also allows a gradual simplification of grain shapes in order to study the importance of the accuracy of the grain representations in the simulation.All this needs to be accomplished at the minimum possible amount of utilized discrete elements in order save on the processing time, which is substantial in a simulation of this type.
Usually, the grains in DEM simulations are represented by simple elements such as spheres or ellipsoids and the irregularity of grain shapes is taken into account through parameters such as rolling friction, that simulate grain interlocking [8].However, estimating the rolling friction solely from grain shapes remains very problematic [9].To directly model the shape of grains, discrete elements can be modelled as polyhedrons -this however proves to be computationally very expensive due to complex contact detection and calculation of contact forces.A very efficient manner of representing arbitrary grain shapes is by clumping together a subset of simpler discrete elements -usually spheres -to approximate the overall grain shape.Hagenmuller et al. [7] used a very straight-forward approach and meshed the grains in the latter manner by simply arranging the spheres on an orthogonal grid.This approach resulted in a large number of discrete elements and very long computational times.In the present article, we present a sphere overlapping grain approximating method that allows us to substantially diminish the number elements needed for grain approximation, while keeping the key geometrical parameters of grains (Figure 1).Although the overlapping sphere technique had already previously been used for grain approximation in several different cases [10,11], our demands on diminishing the number of discrete elements as much a possible while still keeping control over morphological features of grains led us to devise a novel optimized approach.It is based on the concept of medial axis which has been utilized in order to compress the information, needed to reconstruct a binary image [12].
The grain approximating technique is applied by first calculating the medial axis of each grain in order to produce the skeleton of the sample, upon which the approximating spheres will be placed.Placing the center a sphere to every point of the medial axis at this point would result in a perfect reconstruction of the grain image, but would on the other hand contain redundant information due to the discrete nature of the image [12].The medial axis is therefore filtered by removing the spheres which are completely covered by the union of all the other spheres.By starting this filtering process with the smallest spheres and continuing towards the largest sphere, one promotes the use of large spheres for grain reconstruction.
Approximation is carried out via two controlling parameters -the smallest sphere radius R min limits the size of the smallest sphere used for grain approximation and disregards the smaller reconstructing spheres.The delete ratio D deletes all the medial axis points in the distance DR (where R is the radius of the sphere) around the center of an approximated sphere.The spheres are created in the order from the largest to the smallest.Effectively, the parameter D sets the density of the spheres that form the grain surface -the larger the value D, the larger the artificial grain surface roughness, that emerges as a consequence of approximating an originally smooth surface with a set of spheres as shown in the left column of Figure 1.On the other hand, the parameter R min rounds off the sharp edges of the grain and disregards the details on grain surface with a characteristic size smaller that R min as shown in the right column of Figure 1.Hence, one can view the D parameter as the parameter that controls the artificial grain roughness and R min as the parameter of grain detail resolution.

Sensitivity analysis of geometrical features to approximation algorithm parameters
An extensive sensitivity analysis was performed in order to assess the influence of the two grain approximating parameters (D and R min ) on the key geometrical properties of the grains: the volume, the anisotropy and the surface roughness of the grains.
In the present paper, only the study of the volumetric parameter will be presented in details.One is interested in the volumetric error, defined as the volumetric difference between the original and the approximated grains, divided  by the total volume of the original grains.The effect of the two grain approximating parameters on the volumetric error of a subsample of 10 grains is presented in Figure 2, whereas the resulting number of spheres is presented in Figure 3.
As one might expect, decreasing the parameters D and R min will increase the number of spheres, needed for grain approximation (Figure 3), due to a rising density of spheres along the grain surface as well as adding increasingly smaller spheres to represent smaller details on the grain surface.On the other hand, this will also diminish the volumetric error of the approximation (Figure 2).Increasing the parameters D and R min will inevitably produce the opposite effect -diminishing the number of spheres and increasing the volumetric error.

Mechanical behaviour
YADE DEM tool [13] is used to perform the numerical simulation of a classical oedometric compression test on a cubic sample of roundgrained snow with a sidelength of 4,9 mm and initial void ratio of 2.79 (left side of the Figure 4) in order to demonstrate the functioning of the devised micromechanical model.The specimen is approximated with the choice of parameters D = 0.7 and R min = 6, resulting in 19238 discrete elements needed to represent the specimen consisting of 1393 grains.It is confined by a set of rigid walls and compressed with a constant velocity of 10 −2 ms −1 .A cohesive elastic contact law with the brittle regime of cohesion breaking has been used.
Cohesion along the contact surface between the grains is imposed by fist estimating, from the original image, the area of the bond between each pair of grains in contact.Each pixel on the contact surface of a grain is assigned the nearest approximating sphere.Afterwards for each pair of pixels in contact, a cohesive bond is created between the two pertaining approximating spheres.If more than one pair of pixels in contact points to the same pair of approximating spheres, the magnitude of the cohesive bond is multiplied accordingly.The stress-strain curve of the oedometric compression simulation on the snow sample is shown in the Figure 5, together with the corresponding evolution of cohesive grain bond braking (plotted as the ratio between the number of broken bonds and the number of cohesive bonds in the intact sample).A zoom into the primary phase of the curves is provided in the upper left corner of the plot.One can note that the specimen undergoes a short elastic phase in the beginning of the simulation, followed by a peak (ε ≈ 0.001) and subsequent softening 0.001 < ε < 0.01.The amount of broken cohesive bonds remains very low in the elastic phase, rises rapidly at the peak stress, and continues to grow during the softening phase.Further advancement of the mechanical loading results in a multitude of smaller peaks marking the gradual breakage of cohesion, allowing for grain rearrangement and compaction.
The diminishing void ratio results in a gradually increasing resistance to compaction for strains ε > 0.3.The specimen in the final stages of compaction is shown in the right side of Figure 4.
Figure 5.The stress-strain curve of the oedometric simulation on the given snow sample (blue curve) and the evolution of the ratio of broken cohesive bonds between the grains (yellow curve).The corresponding void ratio is marked at the top of the plot.An additional zoom into the primary phase of the simulation is added in the top left corner of the plot.

Conclusion
In the present paper we show development of a micromechanical model, capable of of predicting the mechanical response of a sintered porous material under largestrain loading.The development is based on a pre-existing model, where a novel method of irregular grain representation is added, which results in a substantial diminishment of the number of discrete elements needed for the simulation.An original approach to modelling the cohesion between the grains is also utilized that insures a homogenity of the cohesion along the contact surface of the grains.The parameters of the grain approximation are chosen by performing a thorough sensitivity analysis and studying the effect on the key geometrical features of the approximated grains -the volumetric error of the approximation, the average ratio of the anisotropy of the approximated grains with respect to the anisotropy of the grain images and the amount of the artificial roughness of the approximated grains due to sphere representation of grain surfaces.
The preliminary mechanical evaluation of the model has been conducted in order to demonstrate the functioning of the model on a case of a simple oedometric compression of a snow sample.The results, although preliminary, exhibit all the crucial phases that are to be expected of performing compression on a sintered porous structure: a short elastic phase, followed by a multitude of peaks, marking the gradual breakage of cohesion between the grains.The dimininishig void ratio eventually leads to a rise in mean stresses as the sample gets increasingly compacted.
In our future work, we intend to use this model, in order to study the effect of the microstructure of a weak layer in snow on its mechanical behaviour.Systematic numerical experiments will be conducted on a number of different weak layer sample images in order to determine the dependency between the morphological properties and the peak force and subsequent softening of the weak layer.

Figure 1 .
Figure1.The effect of the two grain approximating parameters.The left column demonstrates the effect of increasing the parameter D, while the right column demonstrates the effect of increasing the parameter R min in terms of approximating grain geometry as well as the number of spheres used for approximation.

Figure 2 .
Figure 2. The effect of the two grain approximating parameters on the volumetric error of approximating the given set of grains.

Figure 3 .
Figure 3.The effect of the two grain approximating parameters on the number of spheres, needed to approximate the given set of grains.

Figure 4 .
Figure 4. Images of the snow sample.The left image shows the sample prior to mechanical loading (ε = 0) while the one on the right shows the sample at the strain of ε = 0.5.