Discrete element simulation of dry snow using the developed analytic bond model

Snow is a heterogenous, hot material which is constituted from ice particles. The bonding behavior of ice particles is an important parameter determining the macroscopic behavior of snow. Discrete Element Method (DEM) is usually used as a tool to model dry snow. The most important input data required into the DEM is bonding behavior of ice particles since ice particles can adhere to form bonds when they brought into contact. This study had two aims: first, an analytical formulation was derived to predict the bond diameter of ice-ice contacts as a function of time, compressive load, and strain rate. Using the previously published data for strain rate of ice, a solution method was developed. The results of bond diameter development with time were compared to experimental data and a good agreement was found. Second, a DEM for dry snow was developed and programmed in MATLAB and the developed bond model was employed in the simulation to study the deposition behavior of snow in a container under gravity acceleration. A specific beam element with implemented damage model was developed in implemented in the simulation using the bond data obtained from the analytical approach. The simulated parameters were macroscopic angle of repose, packing density, and surface conditions as a function of temperature and filling rate. The results showed that discrete element simulations were able to verify the existing published experimental data. Specifically, the simulation results showed that angle of repose of snow decreased rapidly with decreasing the temperature, the surface became very irregular due to the particles rotation and re-arrangement for lower falling speeds of particles, and density increased with depth of deposition. These findings were all matched with experimental observations.

Most particles, including snow, are naturally non-spherical; therefore, it is desired to simulate irregular particles in DEM. Generally, two avenues are paved to treat irregular particles in DEM; one is based on ellipsoids or some polygons and the other route uses agglomerates of bonded spheres/disks. For instance, snow particles can take a plethora of shapes [13], which can be modelled by including bonds between particles [5,6]. If the second approach is chosen, bond modeling is required. To model bonded or cohesive granular materials, traditionally a compact spring model has been used [14][15][16][17][18][19][20]. Despite its wide application, the bending-shear coupling is neglected in most cases. A more sophisticated approach is to employ beam Finite Elements (FEs) to represent the  [21][22][23]. When employing a beam element, it is important to acknowledge that bonds can have a low length to diameter ratio, which makes the Timoshenko beam model a better solution than the Euler beam [22].
Bonds can vary with time, the so called "sintering", which has been studied extensively for metallurgical applications [24]. There has been experimental investigations into the sintering force of ice particles [25,26], and analytical modeling of diffusion sintering [27,28].
The aim of this study is to incorporate a model for short term sintering in DEM and to employ the code to simulate the effects of different parameters including friction coefficient, bonding strength, and deposition rate on the packing configuration of granular materials when short-term sintering is considered.

Contact law for non-bonded particles
The general nomenclature of the contact models for non-bonded and bonded particles are presented in Figure 1. In all the subsequent formulations, the left superscript indicates time or time step, left subscript indicates involved particle(s), right superscript indicates mathematical exponent of the variable, and right subscript indicates abbreviations like " " for rotation or " " for slip. Moreover, italic bold letters represent vectors and non-italic bold letters represent matrices.
In this study, the following damped contact laws are considered for the normal force, , tangential force, , and rolling moment, , respectively [29][30][31] Figure 1. schematic of two particles, and indicating normal penetration, and linear and angular velocity vectors, and , respectively. − is the local normal-tangential coordinate system of the two non-bonded particle system (left), Beam element nomenclature used in this study for bonded particles (right) where, "sign" is the signum function; is the radius of the particle; and are the normal and tangential damping coefficient, respectively; , , and are normal, tangential, and rolling stiffness values, respectively. The expressions for , and can be found in Marshall 2009 [32]; Δ is the relative tangential velocity in direction of the local − coordinate system in Figure 1; is the static friction coefficient; and ̇ are normal penetration and normal penetration speed, respectively. Finally, is a free integration variable.
Relative slip and angular velocities between particles and , at time , are defined as, Δ , Δ , respectively, and implemented as, where, is the angular velocity of particle at time , is the radius of particle , and unless otherwise stated, all units are in SI. is the velocity of particle in direction s at time t. The damping coefficients and have significant effect on the behavior of granular material since they have an important role in dissipating collision energy. Consequently, it is of paramount importance to model a reasonable damping coefficient. A model followed from the work of Tsuji et al. [33] are employed in the current study in which the damping coefficient is considered as a fraction of the critical damping coefficient, 0.9√2 or 0.1√2 , for velocities below and above the sticking velocity, respectively. The density of ice particles is considered a constant value of 917 kg/m 3 since it does not change significantly for the temperature range considered [34].
The time step, , is chosen as, [21,35,36], where is the mass of the smallest particle in the system; is the maximum stiffness in the system, and = 0.025 or 0.2, for bonded and non-bonded cases, respectively.
The friction coefficient of ice, , is modelled as [37], where, Δ is the sliding speed in m/s , = 273.15K, and is the absolute temperature.
The time-dependent sintering force, ( ), can be estimated from the ultimate strength of the bond, , as, where ( ) is the time-dependent bond radius.
The grain strength depends on particle and temperature, as [39], are Boltzmann coefficient, and gas constant, respectively; is absolute temperature; = 3 [39] and R* is again the equivalent radius.
The only remaining parameter that needs to be determined in order to estimate the sintering force is the radius of the contact area, ( ), as depicted in Figure 2.
The initial contact radius, 0 = ( = 0, = 0), when the load is zero [43] can be approximated from JKR theory as, In the case of loaded particles with some external load, , a solution for the viscoelastic case will be possible but it will depend on the loading and unloading history, which is not readily available in applications. As a result, for the initial time-step, the contact radius, = ≠0 ( = 0), of the loaded particles with load are also assumed from JKR, as, Solving for , the initial sintering force when the particles are loaded with a load , can simply be calculated from Eq.
is used [44], where, is applied contact pressure in MPa, and is time in seconds. Assuming that the particle radius, , as depicted in Figure 2, can be considered as a constant during a small time-step Δ , the change in the radius of the contact area can be written as, where is the bond radius at any time , and ̇ is the time derivative of Eq. (15). In practice, this assumption means that negligible melting and diffusion mass transfer from the liquid layer to the bond neck area occurs. Therefore, with / ≪1 for initial sintering stages, 2 − 2 ≈ 2 , and, where ≠0 ( = 0) is the initial radius obtained from Eq. (13). The maximum value of contact pressure, , used in Eq. (17), can be approximated as [45], Finally, the dependence of Young's modulus of the ice to temperature is in all calculations assumed as [46], [GPa] = 9.61(1 + 1.07 × 10 −3 + 1.87 × 10 −6 3 ) −1 (19) where is the temperature in centigrade, and is Young's modulus.
Examples of radius of bond neck area and sintering force versus sintering time under the contact pressures 0.5MPa and 1MPa are presented in Figure 3. It is evident from this figure that with increasing contact pressure, which is equivalent to an increment of the applied load, , the sintering force and the bond radius grow faster. For example, consider a particle of 100 m diameter. For the case with no external loading, the adhesive force using Eq.(14) will be 98 N, which is in agreement with published experimental results [25]. However, for the case where the particle diameter is 6mm, and there is no external loading, the sintering force predicted by Eq. (14) is 3mN, which is an order of magnitude less than that predicted by Szabo and Schneebeli [26]. As depicted in Figure 4(left), the sintering force increases both with time and with particle diameter. The rate of changes of the sintering force is also higher for larger particles, as seen in the figure. The increasing nature of sintering force, also fits well with published experimental results for = 0.5MPa [25]. The sintering force changes with the particle diameter in a nearly linear way if the maximum contact pressure is kept at a constant value, as depicted in Figure 4(right).
where, is the ordinary Timoshenko beam stiffness element, which can be found in any finite element book [47], and is the beam nodal velocity vector of the beam. The Coulomb-Mohr criterion is used to identify the bond failure loads.

Implementation of the method
The algorithm is implemented in MATLAB 2020a ® (Mathworks Inc., Natick, MA). The program have four main functions: 1) generating and filling the particles with location and velocities at top of the container to be thrown into the container, 2) finding possible collisions for each particle at each time step using a sub-cell collision detection method, 3) calculating forces and updating positions using the velocity Verlet method, as presented in Eq.(21) [48], and 4) postprocessing. Positions and velocities are updated as, where, , , and represent displacement, velocity, and acceleration vectors; left superscripts represent the time, and is the marching time from Eq.(6). Similar formulas are used for angular velocities and rotations. Ice particles are generated using a beta-distribution, with the possibility to bond the particles together in the beginning, before inserting them into the container area. At each time step, forces are calculated, and bonds and sticking of particles are checked. Accelerations are obtained simply from forces, and integrated twice to obtain the positions using Eq. (21). The required time for 1000 particles to be simulated for one-million time-steps on a Mac book Pro (CPU 2.2 GHz 6-Core Intel Core i7, 16 GB RAM 2400 MHz DDR4) was 75 minutes.

Numerical experiments
Numerical experiments are conducted by throwing both non-bonded and bonded particles into a container. The effects of friction coefficient, and deposition rate on the packing configuration in either case are investigated. Particles are thrown from the top of the container with an initial random velocity in horizontal direction, and without any initial velocity in vertical direction. Since falling velocity of snow particles is sharply decreasing with their mass or diameter [49] and since the velocity can be in the order of fractions of mm/s for micro-meter sized particles [50], heavy drag forces are introduced for the case of bonded-particles in the current numerical simulation to reduce the effective fall speed. For the non-bonded particles, the simulation time step is in the order of a few micro-seconds and for the bonded particles, it is in the order of a fraction of micro-seconds (≈ 0.3 × 10 −6 s). Five different numerical experiments are detailed and discussed in the following subsections.
A schematic of the numerical experiments of this study are presented in Figure 5, with the critical parameters and variables specified in the figure. In Sec.5.1. , numerical experiments on non-bonded, temperature-independent particles are detailed, while in Sec.5.2. numerical experiment on bonded, temperature-dependent particles are presented. A typical filling process for bonded ice particles, is also shown in Figure 5. Ice particles form in general single particles or 2-, and 3-particle clusters after collisions, like in reported experimental results [50]. After a few particles have filled the bottom, particles start to pile up in a sticky way. In all simulations, the friction coefficient between the ice particles and the container walls (including bottom wall) of the container is set to a constant value of IOP Publishing doi:10.1088/1757-899X/1190/1/012015 8 0.4. A typical plot of particle size distribution, when poly-dispersed particles are considered is presented in Figure 6 ( = 250 + 500β(5,10) m), where β is the beta-distribution [51]. This size distribution is used in all the upcoming simulations considering poly-dispersed particles.

Non-bonded temperature-independent particles
The settings for non-bonded particles simulation are: time average time step, = 1.8 micro-seconds; container size, × ℎ = 20mm × 60mm; particles are released with initial zero velocity; Particles are released in the form of nonbonded linear clusters as depicted in Figure 5 such that the height and horizontal location of the left particle are 54mm and 6.5mm, respectively; cluster length is 0.35 of the width of the container; the nonbonded cluster has a random inclination between 0 (horizontal configuration for the cluster) to 90° (rotated in clockwise direction); Number of particles, = 500; Young's modulus, = 4GPa; Poisson's ratio, = 0.4; Density, = 917Kg/m 3 . Gravitational acceleration and friction coefficient are varied for each experiment and expressed separately.

Mono-dispersed vs. poly-dispersed non-bonded, frictional particles.
In these experiments, gravitational acceleration, = 9.81 / 2 and friction coefficient, = 0.9 are considered and the filling process is according to Figure 5. For experiments with mono-dispersed particles, the IOP Publishing doi:10.1088/1757-899X/1190/1/012015 9 corresponding diameter is set to 500 m and for experiments with poly-dispersed particles the particle diameter varies according to Figure 7. As depicted in Figure 7, deposition of mono-dispersed particles shapes a "preferred"-hexagonal-order, which is the ordered-packing with a preferred highest density arrangement in two dimensional lattices. The compatibility mechanism that allowed the particles to achieve this hexagonal arrangement, is the movement of slip band (shear bands) within the material. Indeed, via geometrically necessary dislocations [52], the material adopted itself for strain gradients during the filling process to take the shape of the container. Voids are also formed as seen in Figure  7(left), which can alter the macroscopic behavior of the material significantly. Adding a perturbation in the size of particles to have nonuniform diameters results in an amorphous configuration as presented in Figure 7(right).

Effects of variation in friction coefficient for non-bonded particles.
Since the static friction coefficient of ice can reach to near unity [37,53], and the related static friction force is the force responsible to allow particles to move and rearrange, an increasing friction coefficient must decrease the packing density, as observed in previous published works [54]. Moreover, it is known that the density of packing increases with the depth of deposition [55].
Here a set of simulation results for non-bonded poly-dispersed particles and various friction coefficients revealed that the density decreases with increasing friction coefficient and increasing with depth, as presented in Figure 8. Gravitational acceleration, = 9.81 / 2 is considered. The left image shows the simulation outcomes while the right image shows the average trends. No temperature dependence is considered on these plots.

Effect of deposition rate on the deposition density for non-bonded, frictional particles.
To attain lower rate in the deposition process, it is possible to suspend the particles in a fluidic medium [56,57]. In this case, due to the movement of particles in the fluid, several forces and moments can be exerted on the particle [32]. For smaller particles, the most important forces are drag forces, and reduced gravity [32].
Here, the effect of deposition rate is studied by reducing the effective gravitational acceleration on the particles. When → 0, particles hit the surface with lower kinetic energy, as presented in Figure  9. It is obvious from this figure that decreasing effective gravity decreases the effective density especially near the surface and increases the irregular pattern depth close to the surface. Gravitational acceleration, = {0.490,0.981,9.81} m/s 2 and friction coefficient = 0.9 are considered while the particles are kept at the constant diameter of 500 m.  Figure 10 that as particles fall, they form new bonds, roll over each other and reconfigure, as previously reported [50]. The settings for this simulation are: gravitational acceleration, = 9,81 / 2 ; poly-dispersed particle size distribution according to Figure 7; average time step, = 0.1 s; container size, × ℎ = 20mm × 60mm; particles are released with initial zero velocity; Particles are released in the form of bonded linear clusters as depicted in Figure 5 such that the height and horizontal location of the middle particle are 54mm and 30mm, respectively; cluster length is 0.35 of the width of the container ; the bonded cluster has a random inclination between 0 (horizontal configuration for the cluster) to 90° (rotated in clockwise direction); Number of particles, = 1000; Poisson's ratio, = 0.4; density, = 917Kg/m 3 ; Friction coefficients, damping values, and Young's modulus, can be extracted based on the information provided in previous sections; drag forces are introduced as = −10 , where and , are velocity and diameter of the particle respectively. Figure 10. experimental results of sticky ice particles accumulation from experimental results [50] (left); simulated results using DEM (right) = −10 Note: For Figure 10, which is copied here from [50], permission will be acquired from Elsevier As a granular material, angles of repose (of snow) is a function of interparticle friction, shape, and bond forces [58]. Therefore, checking the model against an experimentally validated macroscopic model is a reliable way to discern that internal parameters of the model work righteously. A Typical repose configuration of snow at different temperatures are presented in Figure 11. The approximate measured angle of repose from the simulated packing in Figure 11 are compared with the results from Willibald et al [58] and shows good agreement with the experimental results for faceted snow. In Willibald et al. [58] the angle of repose for ice beads are significantly lower which indicates that mono-dispersed particles probably cannot capture the behavior in the best way.
The settings for this simulation are: gravitational acceleration, = 9.81 / 2 ; poly-dispersed particles are considered with an average time step, = 0.1micro-seconds; container size, × ℎ = 60mm × 60mm; particles are released with initial zero velocity; particles are released in the form of bonded linear clusters as depicted in Figure 5 such that the height and horizontal location of the middle particle are 54mm and 30mm, respectively; cluster length is 0.35 of the width of the container; the bonded cluster has a random inclination between 0 (horizontal configuration for the cluster) to 90° (rotating in clockwise direction); Number of particles, = 1000; Poisson's ratio, = 0.4; Density, = 917Kg/m 3 ; friction coefficients, damping values, and Young's modulus, can be extracted based on the information provided in previous section; drag forces are introduced as = −3 , where and , are velocity and diameter of the particle respectively. In contrast to the results in Figure 11, when the container is large enough so there are no vertical wall effects, the relative density decreases with decreasing temperature: ( = −2℃) = 0.68 , while ( = −15℃) = 0.58. Trend-wise, these decrements are is in agreement with the results from experimental works [59]. However, the amount of density decrement is less pronounced than that reported in the literature for natural snow [59], which can be related to the variable natural conditions and wind speeds in nature, which are expected to affect the density significantly. Interestingly, when temperature decreases, angle of repose also decreases, despite the fact that both friction coefficient and bond strength increase. This indicates the importance of bonding in the snow since with decreasing temperature, particles are less likely to stick, according to Eq. (9), and bond growth rate is also slower. In a word, it seems that, the effects of sticking and bond growth dominate the effects of friction and bond strength. To realize the sintering effect in more details, after obtaining the temperature dependent angle of repose for bonded particle, the sintering ability was turned off in the model keeping all other parameters the same, and the simulations was repeated. The results, as depicted in Figure 11(right) show that if the bonding is not considered, the angle of repose is insensitive to the temperature variations.

Conclusions
The discrete element method was used in this study to model the random loose packing behavior of non-bonded and bonded granular materials representing dry snow with different parameters. The effect of variation in friction coefficient and bond strength on random loose packing density was investigated. Bonds were modelled using Timoshenko beam element theory. With parameters described in previous sections, the developed model was verified for a few simple cases where analytic solutions are available. A time-dependent bond model was briefly presented in this study for short-term sintering of granular materials like snow, and the model was validated with previously published experimental results. This study has the following constraints: First of all, only dry snow was considered. Also, fall velocity was considered as a random variable, and, since the focus was on the material model and simpler presentation on the results, simulations were performed in 2D. Moreover, the effect of particles changing size due to the metamorphism was neglected, and the temperature gradient within the accumulated ice was not considered. Furthermore, particles crushing was not considered, and it was assumed that the clusters can break down only to the constitutive particles level. Finally, since we could not find data on the dependency of rolling friction coefficient on the temperature or velocity, it is considered as a constant in this study. If such data are provided, it is of interest to study the effects of variations in rolling friction coefficient on the macroscopic behavior of granular materials like snow. The results showed that in the case of mono-dispersed non-bonded particles, an ordered packing was obtained upon sequential deposition and shear bands or artefacts were formed in the materials to geometrically allow for maximum density and ordered packing. However, adding a small perturbation to the particle size distribution destroyed the ordered packing and an amorphous packing was obtained. An increased friction coefficient decreased the loose packing density of the non-bonded particles. Moreover, it was found that the packing density of the non-bonded particles generally increases with the depth while reducing the deposition velocity decreases the density and generally alters the surface morphology significantly. In packing of bonded particles, deposition velocity also plays an important role. Rolling friction plays a significant role in the packing process and specially in forming the surface morphology. Angle of repose of snow, as a sample bonded material with temperature dependent properties, decreases with decreasing temperature, however becomes asymptotic below -15℃. The angle of repose is insensitive to the temperature, when the bonding ability is not considered. In the case of non-contained packing, where there are no vertical wall effects, the density of snow decreases with the temperature significantly. The results can be useful in the study of random loose packing of snow and any other granular materials and in the study of surface engineering using the deposition method.