# Mechanics of densely packed snow across scales: A constitutive modelling strategy based on the 3-D H-model

Marie Miot, Antoine Wautier, Pit Polfer, Pierre Philippe, Tibor Fülöp, François Nicot

<figure><img src="image_1.png" /></figure>

## To cite this version:

Marie Miot, Antoine Wautier, Pit Polfer, Pierre Philippe, Tibor Fülöp, et al. Mechanics of densely packed snow across scales: A constitutive modelling strategy based on the 3-D H-model. Journal of Glaciology, 2025, 71, pp.e62. 10.1017/jog.2024.112. hal-05457929

HAL Id: hal-05457929

https://hal.inrae.fr/hal-05457929v1

Submitted on 16 Jan 2026

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
Article

**Cite this article:** Miot M, Wautier A, Polfer P, Philippe P, Fülöp T, Nicot F (2025) Mechanics of densely packed snow across scales: A constitutive modelling strategy based on the 3-D H-model. *Journal of Glaciology* **71**, e62, 1–19. https://doi.org/10.1017/jog.2024.112

Received: 8 March 2024
Revised: 17 December 2024
Accepted: 21 December 2024

**Keywords:** constitutive law; multi-scale model; snow

**Corresponding author:** Marie Miot;
Email: miot.marie@yahoo.fr

Marie Miot: Now at University of Calgary,
Calgary, Canada

Mechanics of densely packed snow across
scales: A constitutive modelling strategy
based on the 3-D H-model

Marie Miot¹ ID, Antoine Wautier², Pit Polfer¹, Pierre Philippe², Tibor Fülöp¹ and François Nicot³

¹GIC*L, Goodyear S.A., Colmar-Berg, Luxembourg; ²RECOVER and ISTerre, Aix Marseille Université, Aix-en-Provence, France and ³RECOVER and ISTerre, Université Savoie Mont Blanc, Le Bourget du Lac, France

Abstract

In many technical and geomechanics applications, for example tire and ski design or avalanche prediction, the capability to model the mechanical behaviour of snow is of high importance. To this end, we propose in the present study to extend the 3-D H-model, a multi-scale constitutive law originally developed for granular materials, to densely packed snow. In the model, single ice grains are described by spherical particles bonded by brittle elasto-viscoplastic bridges. Snow is thus described explicitly through its ice skeleton microstructure. As a validation, confined compression test results from the litterature are used to assess the suitability of the model to correctly describe snow behaviour. Multiple parameter studies were conducted to demonstrate the capability of the model to capture the behaviour of different snow types over a significant range of temperatures and loading rates at small deformations. Notably, the initial bond radius emerges as an effective proxy for snow aging under isothermal conditions, with stress levels increasing directly with the initial bond radius. Additionally, low strain rates and elevated temperatures are shown to influence the viscous response of ice bonds, their failure rates and the overall stress within the snow material.

**1. Introduction**

As a granular material constituted of an interconnected network of ice grains, snow has a complex mechanical behaviour influenced both by the material properties of ice and by its intrinsic microstructure. As a homogeneous material, ice can show both brittle and ductile behaviour, depending on the external conditions as strain rate for instance (Schulson, 1990). It also exhibits temperature-dependant creep behaviour (Barnes and others, 1971). From a microscopic perspective, snow exhibits diverse microstructures initially dictated by the meteorological conditions at snowfall but evolving over time under the influence of the environmental constraints. This so-called snow metamorphism is mostly driven by the temperature and its gradient within the snow layer. Under specific climatic conditions, snow metamorphism leads to the transformation of snowflakes into bonded spherical ice grains (Fierz, 2009; Calonne and others, 2014) which normally goes along with a snow densification (Kojima, 1967).

Modelling the mechanical behaviour of snow is of great interest in various fields, from avalanche prediction and mitigation (Nicot, 2004; Gaume and others, 2015), winter sports equipment (Fauve and others, 2008), to snow tire design (Shoop and others, 2006; Lee, 2013; Hjort and others, 2017). Two kinds of approaches are widely adopted for snow modelling, namely continuum-based and particle-based methods. Particle-based methods have the advantage of explicitly accounting for the characteristics of snow microstructure. They are based on an individual representation of grains, mostly through the discrete element method (DEM) (Cundall and Strack, 1979; Radjai and others, 1998; Luding, 2008; Smilauer, 2010). DEM has already been used for snow by modelling the ice grains as discrete particles connected by ice bonds (Johnson and Hopkins, 2005; Mede and others, 2018; Mulak and Gaume, 2019; Willibald and others, 2019; Bobillier and others, 2021; Kabore and others, 2021; Peters and others, 2021) or to represent snow elements (Gaume and others, 2015). In return for taking effective account of the effect of the microstructure, DEM is computationally costly, which limits its ability to model large scale systems. As an alternative to discrete numerical methods, continuum-based approaches of snow include the finite element method (Meschke and others, 1996; Podolskiy and others, 2013; Fourteau and others, 2024), the finite difference method (Dent and Lang, 1980), the material point method (Gaume and others, 2018, 2019; Trottet, 2022) and the smooth particles hydrodynamics (El-Sayegh and El-Gindy, 2019) among others. In these methods, snow

© The Author(s), 2025. Published by
Cambridge University Press on behalf of
International Glaciological Society. This is an
Open Access article, distributed under the
terms of the Creative Commons Attribution
licence (http://creativecommons.org/licenses/
by/4.0), which permits unrestricted re-use,
distribution and reproduction, provided the
original article is properly cited.
is modelled as a continuous medium and its material behaviour is accounted for by a constitutive relation. For example, a Drucker-Prager (Haehnel and Shoop, 2004; Lee, 2009; Gaume and others, 2011; Blatny and others, 2023) or a modified Cam-Clay (Gaume and others, 2018; Guillet and others, 2023) constitutive material laws are widely used in the field of snow avalanche modelling. Such approaches are usually significantly less time-consuming than DEM, meaning that they can address large scale systems. However, the constitutive relations implemented are generally phenomenological and therefore do not explicitly account for the snow microstructure, which makes them harder to calibrate. However, recent works have highlighted the importance of incorporating snow microstructure into constitutive models to more accurately capture its mechanical behaviour (Blatny and others, 2021, 2023).

In this paper, we propose a continuum approach based on a constitutive material model that incorporates explicitly the discrete microstructure of snow, thus combining the advantages of both modelling strategies. To this end, the 3-D H-model (Nicot and Darve, 2011; Xiong and others, 2017) is extended to snow, considered as an assembly of spherical ice particles connected by ice bonds. This model was originally introduced for soil modelling in dry (Wautier, 2021) and partially saturated conditions (Xiong, 2021). The ability of the H-model to be used at large scale has been demonstrated by Xiong and others (2019), where the 3-D H-model is used to model geotechnical problems. In this approach, a granular material is described as a distribution of mesostructures formed by a few grains. The collection of these mesostructures (denoted H-cells in the following) in various configurations provides a statistical description of the microstructure. From the individual mechanical response of each mesostructure, the material constitutive relation can be formulated in a continuum framework by statistical homogenization. In the hereafter presented extension of the 3-D H-model, the particles forming the H-cells represent ice grains with an interparticle contact law specifically suited to ice. The latter was developed to capture the ice bond evolution as well as the ice grain contact and is partly based on a previous work (Kabore and others, 2021). As many engineering applications (vehicles moving on snow, winter sports, ...) imply short-term time scales and fast loading, only the brittle behaviour of ice will be considered in the following. Since density is one of the most important parameters in snow behaviour, special attention is brought to the relationship between the geometry of the mesostructures in the H-model and the macroscopic density. The influence of the aging of snow before the mechanical loading and its impact on the calibration of the numerical parameter are also investigated. Finally, the capacity of the model to account for varying conditions such as the temperature and the loading rate is demonstrated.

## 2. Multi-scale model for snow

The H-model was initially developed for granular soils, first in two dimensions (Nicot and Darve, 2011; Veylon and others, 2018) before being extended to three dimensions (Xiong and others, 2017). It has been progressively extended to capture additional physical phenomena, such as the existence of capillary forces in unsaturated soils (Wautier, 2021; Xiong, 2021), or erosion by suffusion (Ma and others, 2022). Here, further extension of the 3-D H-model to snow modelling relies on implementing a specific customized contact law for ice bonds, restricting the scope of the study to fast loading conditions (strain rate above $10^{-4} s^{-1}$), so that the ice behaviour can be assumed to be exclusively brittle (Schulson, 1990).

<figure><img src="image_3.png" /><figcaption>**Figure 1.** 3-D bi-hexagonal mesostructure in the global reference frame with the definition of the Euler angles.</figcaption></figure>

### 2.1. A multi-scale constitutive model for granular materials: the 3-D H-model

Contrary to phenomenological approaches, multi-scale models introduce a statistical homogenization process, which consists in inferring the mechanical behaviour of a given material from the constitutive properties taking place at a smaller microscopic scale (Wautier, 2021). More precisely, the material constitutive law is obtained from the weighted averaged behaviours of several mesostructures, which are each constituted by a set of grains. First conceptualized by Satake (1992) and later explored by Kruyt and Rothenburg (1996), grain loops are an example of such mesostructures in a 2-D granular assembly. In principle, these loops enable a geometrical description through a particle graph paving the space and have been shown to be a powerful tool to explain the behaviour of granular materials based on local physics. Namely, it has been observed that grain loops containing more than six grains within a granular material can serve as a meaningful descriptor for the material scale behaviour (Zhu and others, 2016, 2016a, 2016b; Liu and others, 2018). These results support the physical suitability of the H-model in which the microstructure is described as a collection of either hexagonal grain assemblies in the 2-D H-model (Nicot and Darve, 2011) or bi-hexagonal grain assemblies in the 3-D H-model (Xiong and others, 2017).

The tridimensional bi-hexagonal mesostructure is formed of two perpendicular hexagons as represented in Fig. 1. Each H-cell has an orientation in space, characterized by the three Euler angles $\phi$, $\theta$ and $\psi$ with respect to the global reference frame. The distribution of this orientation $\omega(\theta, \phi, \psi)$ describes the initial anisotropy of the microstructure. The two hexagons of the mesostructure are supposed to remain in orthogonal planes over
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Schematic representation of the two hexagonal grain configurations constituting the 3-D H-cell showing: (a, d) geometry of the hexagons, (b, e) external forces on the hexagons and (c, f) equilibrium of forces on the grains. The dashed-dotted green lines represent the symmetry axis in the considered plane.</figcaption></figure>

any loading. The planar geometry of the two hexagons is shown in Fig. 2. By adopting symmetry assumptions in the cell geometry, each hexagon can be parameterized by only two intergranular distances ($d_1$ and $d_2$ in hexagon A, and $d_3$ and $d_4$ in hexagon B) and an opening angle ($\alpha_1$ in hexagon A and $\alpha_2$ in hexagon B). The symmetries in the cell allow to consider only four different contact points, between grains 1–2, 2–3, 1–7 and 7–8. To maintain the overall symmetry assumption in the cell geometry, equality is also imposed between the various forces within the cell: $F_i$ and $G_i$, representing the normal and tangential external forces in the three directions of the space, and $N_i$ and $T_i$, representing the normal and tangential contact forces at the four different contact points between the grains, as described in Fig. 2. The external

vertical force on grain 1 can be decomposed into two terms $F_1^a$ and $F_1^b$, which respectively insure the equilibrium of hexagons A and B.

Depending on the loading of the cells, some contacts between grains can appear or be lost. Appendices C and D describe these potential cases in detail.

The overall homogenization scheme used for the H-model, as summarized in Fig. 3, allows connecting continuum/macroscale and particle/microscale descriptions of the material. An incremental constitutive relation linking the global incremental strain tensor $\delta\varepsilon$ to the global incremental stress $\delta\sigma$ is thus implicitly derived. In the following, we briefly review the main steps of this homogenization scheme.
<figure><img src="image_5.png" /><figcaption>Figure 3. Schematic representation of overall sequence of the 3-D H-model methodology.</figcaption></figure>

The model assumes strain homogeneity such that the same macroscopic strain is applied in all H-cells. Consequently, the local incremental strain tensor $\delta\tilde{\epsilon}$ in each H-cell is expressed in the local frame as

$$ \delta\tilde{\epsilon} = \underline{\underline{P}} \delta\epsilon = \underline{\underline{P}}^{-1} \quad (1) $$

with $\underline{\underline{P}}$ being the rotation matrix between the global and local frames:

$$ \underline{\underline{P}} = \begin{pmatrix} 1 & 0 & 0 \\ 0 & \cos \psi & -\sin \psi \\ 0 & \sin \psi & \cos \psi \end{pmatrix} \begin{pmatrix} \cos \varphi & 0 & -\sin \varphi \\ 0 & 1 & 0 \\ \sin \varphi & 0 & \cos \varphi \end{pmatrix} \times \begin{pmatrix} 1 & 0 & 0 \\ 0 & \cos \theta & -\sin \theta \\ 0 & \sin \theta & \cos \theta \end{pmatrix}. \quad (2) $$

The changes in the H-cell lengths $\delta L_i$ can be deduced from the diagonal coefficients of the local incremental strain tensor and the external length of the H-cell $L_i$, as defined in Fig. 2:

$$ \delta L_i = \delta \tilde{\epsilon}_{ii} L_i. \quad (3) $$

The evolution of the H-cell geometry depends directly on the variations of the H-cell lengths and on the force balance applied to the mesostructure. The geometrical compatibility equations read:

$$ \left\{ \begin{array}{l} L_1 = l_1 + 2r_g = d_2 + 2d_1 \cos \alpha_1 + 2r_g = d_4 + 2d_3 \cos \alpha_2 + 2r_g \\ L_2 = l_2 + 2r_g = 2d_1 \sin \alpha_1 + 2r_g \\ L_3 = l_3 + 2r_g = 2d_3 \sin \alpha_2 + 2r_g \end{array} \right. \quad (4) $$

with $L_i$ the lengths of the smallest peripheral parallelepiped volume, $l_i$ the lengths of the hexagons defined by the centres of the grains and $r_g$ the grain radius.

Differential variation of Eqn (4) provides three equations relating the variations of the H-cell lengths $\delta L_i$ to the variations of both the intergranular distances $\delta d_i$ and the opening angles $\delta \alpha_i$:

$$ \left\{ \begin{aligned} \delta L_2 &= \delta d_2 + 2\delta d_1 \cos \alpha_1 + d_1 \sin \alpha_1 \delta \alpha_1 \\ &= \delta d_4 + 2\delta d_3 \cos \alpha_2 + d_3 \sin \alpha_2 \delta \alpha_2 \\ \delta L_2 &= 2\delta d_1 \sin \alpha_1 + d_1 \cos \alpha_1 \delta \alpha_1 \\ \delta L_3 &= 2\delta d_3 \sin \alpha_2 + d_3 \cos \alpha_2 \delta \alpha_2 \end{aligned} \right. . \quad (5) $$

Considering the balance of forces and momentum for the grains 2 and 7 provides:

$$ \left\{ \begin{aligned} G_2 &= T_1 \\ G_3 &= T_3 \\ N_2 &= N_1 \cos \alpha_1 + T_1 \sin \alpha_1 + G_2 = N_1 \cos \alpha_1 + T_1(1 + \sin \alpha_1) \\ N_4 &= N_3 \cos \alpha_2 + T_3 \sin \alpha_2 + G_4 = N_3 \cos \alpha_2 + T_3(1 + \sin \alpha_2) \end{aligned} \right. . \quad (6) $$

Differentiating these last two equilibrium equations yields:

$$ \left\{ \begin{aligned} \delta N_1 \cos \alpha_1 - N_1 \sin \alpha_1 \delta \alpha_1 + \delta T_1 (1 + \sin \alpha_1) + T_1 \cos \alpha_1 \delta \alpha_1 \\ - \delta N_2 = 0 \\ \delta N_3 \cos \alpha_2 - N_3 \sin \alpha_2 \delta \alpha_2 + \delta T_3 (1 + \sin \alpha_2) + T_3 \cos \alpha_2 \delta \alpha_2 \\ - \delta N_4 = 0 \end{aligned} \right. . \quad (7) $$

The contact forces are related to the H-cell geometrical parameters through a contact law. In the general case of a visco-elastic contact, a generic form of the contact law at a given contact i reads:

$$ \left\{ \begin{aligned} N_i &= a_{n,i} (d_i - 2r_g) + c_{n,i} (t - t_{\text{ref}}) \\ T_i &= a_{t,i} (d_i - 2r_g) + b_{t,i} (\alpha_j - \alpha_{j,\text{ref}}) + c_{t,i} (t - t_{\text{ref}}) \end{aligned} \right. \quad (8) $$

with $a_{n,i}$, $c_{n,i}$, $a_{t,i}$, $b_{t,i}$ and $c_{t,i}$ being independent of the geometry. The $a_{n,i}$ and the $a_{t,i}$ coefficients are stiffness terms expressing the evolution of the intergranular distance in terms of the
incremental contact forces (normal and tangential respectively). The coefficients $b_{t,i}$ link the tangential contact force to the evolution of the opening angle. The coefficients $c_{n,i}$ and $c_{t,i}$ are viscosity terms, which depend on the contact loading. Note that $\alpha_{i,ref}$ denotes the opening angle at the beginning of the loading at time $t_{ref}$. The contact law specifically developed for snow will be detailed in the next section.

Differentiation of Eqn (8) provides a generic incremental form of the contact law:

$$ \left\{ \begin{array}{l} \delta N_i = a_{n,i} \delta d_i + c_{n,i} \delta t \\ \delta T_i = a_{t,i} \delta d_i + b_{t,i} \delta \alpha_j + c_{t,i} \delta t \end{array} \right. . \quad (9) $$

Combining Eqns (5), (7) and (9) results in two sets of equations that provide the changes in the intergranular distances and in the opening angles as a function of the changes in the H-cell lengths:

$$ \begin{pmatrix} 2 \cos \alpha_1 & 1 & -2d_1 \sin \alpha_1 \\ 2 \sin \alpha_1 & 0 & 2d_1 \cos \alpha_1 \\ A_1 & B_1 & C_1 \end{pmatrix} \begin{pmatrix} \delta d_1 \\ \delta d_2 \\ \delta \alpha_1 \end{pmatrix} = \begin{pmatrix} \delta L_1 \\ \delta L_2 \\ \lambda_1 \end{pmatrix} \quad (10) $$

and

$$ \begin{pmatrix} 2 \cos \alpha_2 & 1 & -2d_3 \sin \alpha_2 \\ 2 \sin \alpha_2 & 0 & 2d_3 \cos \alpha_2 \\ A_2 & B_2 & C_2 \end{pmatrix} \begin{pmatrix} \delta d_3 \\ \delta d_4 \\ \delta \alpha_2 \end{pmatrix} = \begin{pmatrix} \delta L_1 \\ \delta L_3 \\ \lambda_2 \end{pmatrix} \quad (11) $$

with

$$ \left\{ \begin{array}{l} A_1 = a_{n,1} \cos \alpha_1 + (1 + \sin \alpha_1) a_{t,1} \\ B_1 = -a_{n,2} \\ C_1 = -N_1 \sin \alpha_1 + T_1 \cos \alpha_1 + (1 + \sin \alpha_1) b_{t,1} \\ \lambda_1 = -c_{n,1} \delta t \cos \alpha_1 - c_{t,1} \delta t + c_{n,2} \delta t \end{array} \right. \quad (12) $$

and

$$ \left\{ \begin{array}{l} A_2 = a_{n,3} \cos \alpha_2 + (1 + \sin \alpha_2) a_{t,3} \\ B_2 = -a_{n,4} \\ C_2 = -N_3 \sin \alpha_2 + T_3 \cos \alpha_2 + (1 + \sin \alpha_2) b_{t,2} \\ \lambda_2 = -c_{n,3} \delta t \cos \alpha_2 - c_{t,3} \delta t + c_{n,4} \delta t \end{array} \right. . \quad (13) $$

The incremental contact forces can be obtained from the incremental intergranular distances and opening angles according to the incremental contact law (9).

Continuing with the continuum formulation, the Love-Weber formula (Love, 1892) expresses the local stresses $\tilde{\sigma}$ at the scale of the H-cell from the updated contact forces $N_i(t) = N_i(t - \delta t) + \delta N_i$ and $T_i(t) = T_i(t - \delta t) + \delta T_i$

$$ \left\{ \begin{array}{l} \tilde{\sigma}_{11} = \frac{1}{V_{\text{meso}}} (4N_1 d_1 \cos^2 \alpha_1 + 4T_1 d_1 \cos \alpha_1 \sin \alpha_1 + 2N_2 d_2 \\ \qquad + 4N_3 d_3 \cos^2 \alpha_2 + 4T_3 d_3 \cos \alpha_2 \sin \alpha_2 + 2N_4 d_4) \\ \tilde{\sigma}_{22} = \frac{1}{V_{\text{meso}}} (4N_1 d_1 \sin^2 \alpha_1 - 4T_1 d_1 \cos \alpha_1 \sin \alpha_1) \\ \tilde{\sigma}_{33} = \frac{1}{V_{\text{meso}}} (4N_3 d_3 \sin^2 \alpha_2 - 4T_3 d_3 \cos \alpha_2 \sin \alpha_2) \\ \tilde{\sigma}_{ij} = 0 \text{ if } i \neq j \end{array} \right. \quad (14) $$

where $V_{\text{meso}}$ is the external volume of the mesostructure (Wautier, 2021). Finally, a statistical homogenization of the stresses over the complete H-cell collection provides the stress tensor $\underline{\sigma}$ in the global frame:

$$ \underline{\sigma} = \frac{1}{V_{\text{macro}}} \iiint \omega(\theta, \varphi, \psi) V_{\text{meso}}(\theta, \varphi, \psi) \underline{P}^{-1} \underline{\sigma}(\theta, \varphi, \psi) \times \underline{P} \sin \varphi \, d\varphi \, d\theta \, d\psi \quad (15) $$

where $\omega(\theta, \varphi, \psi)$ is the probability density function of the H-cell distribution, which satisfies:

$$ \iiint \omega(\theta, \varphi, \psi) \sin \varphi \, d\varphi \, d\theta \, d\psi = 1 \quad (16) $$

and $V_{\text{macro}}$ is the total volume of the distribution of H-cells given by:

$$ V_{\text{macro}} = \iiint \omega(\theta, \varphi, \psi) V_{\text{meso}}(\theta, \varphi, \psi) \sin \varphi \, d\varphi \, d\theta \, d\psi. \quad (17) $$

Eventually, we obtain a multi-scale constitutive model that relates the total stress to the global strain under an incremental formalism (Fig. 3). The hypothesis of non-correlated mesostructures inferred by the statistical homogenization, and the hypothesis on the geometry of the mesostructure could lead to an error, but previous works on the H-model (Nicot and Darve, 2011; Xiong and others, 2017; Veylon and others, 2018; Wautier, 2021; Xiong, 2021; Deng, 2022; Liu and others, 2024) have proven that this error is not significant.

## 2.2. Contact law between ice grains

In the extended version of the 3-D H-model proposed hereafter, the snow constitutive features are accounted for by introducing a specific contact law between the ice grains of each hexagon of a H-cell. In the context of the snow-H-model, the history of a contact is supposed to be the following:

(i) The ice grains have undergone metamorphism processes which progressively round their shape (Fierz, 2009). Ice grains are thus assumed to be spherical.

(ii) Thermodynamic sintering has created viscous ice bonds at all contacts between ice grains (Kabore and others, 2021) (Peters and others, 2021). Each ice bond is assumed to be incompressible (Lipovsky, 2022), meaning that its volume originates only from the prior sintering conditions, such as temperature and sintering time (Szabo and Schneebeli, 2007).

(iii) Bond creation processes are supposed to be slow compared to the loading rates considered in the targeted applications ($\dot{\varepsilon} \ge 10^{-4} \text{s}^{-1}$). Consequently, if an ice bond fails during the stress loading, it has no time to reform, and the contact becomes cohesionless until the end of the loading.

Based on these assumptions, the contact law to be developed in the snow-H-model must account for both bonded and unbonded contacts between spherical ice grains.

Grains in contact are assumed to be initially bonded by a cylindrical ice bridge resulting from prior snow metamorphism (Fig. 4). The bridge is characterized by its radius $r_b$ and its length $L_b$, with $r_{b,0}$ and $L_{b,0}$ their respective values at the initial time step. Note that the initial bond geometry before loading created by prior thermodynamic sintering can be changed through its initial radius, enabling the modelling of different snow types (Szabo and Schneebeli, 2007). During mechanical loading, the strain is assumed to be localized into the cylindrical ice bonds while the spherical ice grains are assumed perfectly rigid so that the strain is fully concentrated within the bond.

The deformable bond between two ice grains is modelled by a Maxwell elasto-viscoplastic beam on the basis of the work of Kabore and others (2021) and Peters and others (2021). The main discrepancy with this model lies in the absence of grain rotation in the H-model which leads to the lack of torsion and bending in the ice beam. The creep behaviour of ice, which can be dominant during loading, is accounted for by introducing viscoplasticity in the bond.
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Schematic representation of a bonded contact: (a) initially after snow metamorphism and before loading; (b) during the loading. The yellow hatched zones represent the deformable ice bond, and the blue dotted zones represent the non-deformable parts of the ice grains.</figcaption></figure>

The total normal bond strain $\epsilon_{b,n}^{tot}$ can be split into an elastic part $\epsilon_{b,n}^{e}$ and a viscoplastic part $\epsilon_{b,n}^{vp}$:

$$\epsilon_{b,n}^{tot} = \epsilon_{b,n}^{vp} + \epsilon_{b,n}^{e} \qquad (18)$$

The viscoplastic deformation accounts for the well-known creep behaviour of ice since the creep rate $\dot{\epsilon}_{b,n}^{vp}$ in the bond direction is related to the normal stress $\sigma_{b,n}$ by Glen's law (Barnes and others, 1971):

$$\dot{\epsilon}_{b,n}^{vp} = \frac{\delta \epsilon_{b,n}^{vp}}{\delta t} = A(T_{\text{emp}}) \sigma_{b,n}^{m} \qquad (19)$$

where $A(T_{\text{emp}})$ is a temperature dependent coefficient, and $m$ is an exponent generally close to 3 (Barnes and others, 1971).

Finally, the elastic normal deformation is simply proportional to the normal stress. Due to the adhesion of the bond on the ice grains (which prevents the bond from expanding laterally around the contact points), using Young's modulus $E$ would underestimatethe stiffness of the bond/grains assembly. Thus, the oedometer coefficient of ice $E_{\text{oedo}}$ is used instead to better reflect the lateral geometrical constraints:

$$\sigma_{b,n} = \frac{E(1-\nu)}{(1+\nu)(1-2\nu)} \epsilon_{b,n}^{e} = E_{\text{oedo}} \epsilon_{b,n}^{e} \qquad (20)$$

with $\nu$ the Poisson's ratio of ice.

The normal strain in the bond can be expressed as a function of the bond length, thanks to the strain localization assumption:

$$\epsilon_{b,n}^{tot} = -\frac{L_b - L_{b,0}}{L_{b,0}} \qquad (21)$$

with $L_b$ and $L_{b,0}$ being the bond length, respectively, at the current time and at the beginning of the loading. It should be noted that compression and contraction are counted positive.

The normal stress can be related to the normal contact force $N^b$ as follows:

$$\sigma_{b,n} = \frac{N^b}{S_b} = \frac{N^b}{\pi r_b^2} \qquad (22)$$

where $r_b$ is the radius of the cylindrical bond and $S_b$ its section.

As the bond is the only deformable element in the bonded ice grain model, the incremental bond length $\delta L_b$ is directly equal to the incremental intergranular distance $\delta d$:

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Different modes of bond failure (tensile, compression and shear). Note that potential additional failure modes by bending and torsion are not accounted for as rotation of the grains is not possible in the H-model.</figcaption></figure>

$$\delta L_b = \delta d \qquad (23)$$

Eventually, the incompressibility condition, together with the volume conservation of the ice bond, read:

$$\pi r_b^2 L_b = \pi r_{b,0}^2 L_{b,0} \qquad (24)$$

By combining Eqns (18)–(24), the contact law relating the incremental normal contact force $\delta N^b$ to the incremental intergranular distance $\delta d$ can be obtained:

$$\delta N^b = -\pi r_b^2 \left( \frac{\sigma_{b,n}}{L_b} + \frac{E_{\text{oedo}}}{L_{b,0}} \right) \delta d - \pi r_b^2 E_{\text{oedo}} A(T_{\text{emp}}) \sigma_{b,n}^m \delta t \qquad (25)$$

Likewise, the incremental tangential contact force can be expressed as a function of the incremental intergranular distance and of the variation of the bond direction, characterized by $\delta \alpha$:

$$\delta T^b = -\frac{T \delta d}{L_b} + \pi r_b^2 G \frac{L_b \delta \alpha}{L_{b,0}} - \pi r_b^2 G A(T_{\text{emp}}) \sigma_{b,t}^m \delta t \qquad (26)$$

with $G = \frac{E}{2(1+\nu)}$ the shear modulus of ice. The way to obtain this expression using Glen's law is detailed in Appendix A.
During the loading, stresses increase in the bond, which can lead to its failure. In our approach, only brittle failure is accounted for as the strain rate in the targeted engineering applications is supposed to be high enough ($\dot{\epsilon} \ge 10^{-4} \text{s}^{-1}$). Brittle bond failure occurs as soon as one of the following conditions is met:

$$ \left\{ \begin{array}{l} \sigma_{\text{b,n}} > \sigma_{\text{bcs}} \\ \sigma_{\text{b,n}} \leftarrow -\sigma_{\text{ts}} \\ |\sigma_{\text{b,t}}| > \sigma_{\text{ss}} \end{array} \right. \qquad (27) $$

where $\sigma_{\text{b,n}}$ is the normal stress in the bond, which is positive in compression, $\sigma_{\text{b,t}}$ is the shear stress in the bond. $\sigma_{\text{bcs}}$, $\sigma_{\text{ts}}$ and $\sigma_{\text{ss}}$ are respectively the compression, tensile and shear strength of ice (see Figure 5). The expression of the brittle strengths of the ice bonds is given by Michael (2014) and Kabore and others (2021):

$$ \left\{ \begin{array}{l} \sigma_{\text{bcs}} = 10 \text{ MPa} - \frac{3}{8} T_{\text{emp}} (^{\circ}\text{C}) \\ \sigma_{\text{ts}} = 1.1 \text{ MPa} \\ \sigma_{\text{ss}} = \frac{1}{3} \sigma_{\text{bcs}} \end{array} \right. . \qquad (28) $$

When bond failure occurs, the contact law between the ice grains is reduced to a simple elasto-plastic law:

$$ \left\{ \begin{array}{l} \delta N_i = \begin{cases} -k_n \delta d_i & \text{if } d_i \le 2r_g \\ 0 & \text{else} \end{cases} \\ \\ \delta T_i = \begin{cases} k_t \delta d_i \delta \alpha_j & \text{if } |T_i + k_t d_i \delta \alpha_j| < |N_i + \delta N_i| \\ \text{sign}(T_i + k_t d_i \delta \alpha_j) (N_i + \delta N_i) \tan \varphi_g - T_i & \text{else} \end{cases} \end{array} \right. . \qquad (29) $$

In the absence of ice bond, the normal stiffness $k_{n,i}$ between ice grains at contact *i* evolves with grain interpenetration (see justification in Appendix B):

$$ k_{n,i} = E_{\text{ice}} \frac{\pi}{2} \left( r_g - \frac{d_i^2}{4r_g} \right) . \qquad (30) $$

## 2.3. Extension of the 3-D H-model to snow

The relationship between the incremental contact forces and the changes in the local kinematics (including the intergranular distances between two ice grains and the contact direction) has been previously established in Section 2.2, both for a bonded contact (Eqns (25) and (26)) and after failure of the ice bond (Eqn (29)). The incremental normal contact force in Eqns (25) and (29) and the incremental tangential force in Eqns (26) and (29) can be written generically as in Eqn (9), with the corresponding coefficients:

$$ \left\{ \begin{array}{l} \alpha_{n,i} = -\pi r_b^2 \left( \frac{\sigma_{b,n}}{L_b} + \frac{E_{0ed0}}{L_{b,0}} \right) I_{b,i} - E_{\text{ice}} \frac{\pi}{2} \left( r_g - \frac{d_i^2}{4r_g} \right) (1 - I_{b,i}) \\ c_{n,i} = -\pi r_b^2 E_{0ed0} A (T_{\text{emp}}) \sigma_{b,n}^m \delta t I_{b,i} \\ a_{t,i} = \frac{-T}{L_b} I_{b,i} + \text{sign}(T_i + k_t d_i \delta \alpha_j) E_{\text{ice}} \frac{\pi}{2} \left( r_g - \frac{d_i^2}{4r_g} \right) (1 - I_{b,i}) \\ \times (1 - I_{e,i}) \\ b_{t,i} = \pi r_b^2 G \left( \frac{L_b}{L_{b,0}} \right) I_{b,i} + E_{\text{ice}} \frac{\pi}{2} \left( r_g - \frac{d_i^2}{4r_g} \right) d_i (1 - I_{b,i}) I_{e,i} \\ c_{t,i} = -\pi r_b^2 G A (T_{\text{emp}}) \sigma_{b,t}^m \delta t I_{b,i} + (\text{sign}(T_i + k_t d_i \delta \alpha_j) N_i \\ \times \tan \varphi_g - T_i) (1 - I_{b,i}) (1 - I_{e,i}) \end{array} \right. \qquad (31) $$

where *i* stands for the number of the considered contact, as defined in Fig. 2. $I_{b,i}$ is a bond indicator equal to one when the ice bond

<table><caption>**Table 1.** Experimental (left) and numerical parameters (right) for confined compression tests on snow (Abele and Gow, 1976)</caption><thead><tr><th>Experimental parameters</th><th></th><th>Numerical parameters</th><th></th></tr></thead><tbody><tr><td>Temperature (°C)</td><td>-34; -4; -1</td><td>Temperature (°C)</td><td>-34, -1</td></tr><tr><td>Strain rate (s<sup>-1</sup>)</td><td>5.26; 7.84; 16.0</td><td>Strain rate (s<sup>-1</sup>)</td><td>5.26; 7.84; 16.0</td></tr><tr><td>Initial snow density (kg m<sup>-3</sup>)</td><td>[280; 760]</td><td>Initial opening angle (°)</td><td>[45; 90]</td></tr><tr><td>Snow age (days)</td><td>0.1; 3; 7</td><td>Initial intergranular distance</td><td>≥ 2r<sub>g</sub></td></tr><tr><td>Sample volume (cm<sup>3</sup>)</td><td>[300; 6800]</td><td>Initial bond radius</td><td>[0; r<sub>g</sub>]</td></tr><tr><td></td><td></td><td>Number of H-cell directions</td><td>500</td></tr></tbody></table>

at contact *i* exists, and zero otherwise. $I_{e,i}$ is a similar indicator which is equal to one in a nonbonded elastic contact, and zero otherwise. The ice contact law defined by these coefficients can be implemented readily within the standard computation scheme of the H-model. Some additional improvements of the model in relation with the creation and the loss of contacts in the H-cell are described in Appendix C and D.

## 3. Homogeneous compression test simulations

In this section, the snow-H-model is validated at the representative elementary volume (REV) scale against experimental tests by (Abele and Gow, 1976) conducted for different snow types. The numerical parameters are tuned to match the experimental stress/strain curves as closely as possible. This calibration procedure enables to identify the snow density range over which the model can capture the mechanical response of snow.

### 3.1. Confined compression test

In Abele and Gow (1976), several confined compression tests were performed on different snow samples. In practice, natural snow was collected in cylindrical containers and compacted at temperatures ranging between -7°C and -1°C. The snow samples were then stored for 0.1-7 days before being tested under oedometer conditions at a constant temperature between -34°C and -1°C. The storage duration prior to testing is referred to as the snow age in the following. During the compression test, a vertical velocity of 40 cm/s was applied to the top of the sample, while lateral deformations were kept constant. The experimental parameters are summarized in the first two columns of Table 1. The initial snow density of each sample was determined from the original volume and after measurement of the final sample weight.

### 3.2. Calibration against experimental results

The extension of the 3-D H-model to snow is implemented in a stand-alone Python code which is used to run REV scale tests under strain-controlled conditions. The compression tests under lateral confining from Abele and Gow (1976) are reproduced with this code.

Some of the experimental parameters of the experiments such as the snow temperature and the strain rate can be directly used in the simulations. However, the initial snow density and the snow age are not explicit input parameters for the 3-D H-model. The experimental snow age is only indirectly related to the initial bond radius $r_{b,0}$ (Szabo and Schneebeli, 2007) and the initial density derives from the initial geometry of the H-cells. The latter is determined by the initial opening angle $\alpha_{1,0} = \alpha_{2,0} = \alpha_0$ which can be varied to modify the initial density. Unlike the original H-model for
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Initial configurations of a hexagonal mesostructure: (a) 2-D representation of a ten-grain H-cell with $d_{2,0} = d_{1,0} = 2r_g$ and initial ice bonds bond1 and bond2; (b) 2-D representation of two separated five-grain mesostructure with $d_{1,0} = 2r_g$, $d_{2,0} > 2r_g$ and initial bonds bond1.</figcaption></figure>

sand, some grains may not necessarily be in contact in the snow-H-model (see Fig. 6), to account for large porosity values in snow. A new parameter is then introduced: the initial inter-granular distance $d_{2,0}$ (see Fig. 6 where $d_{2,0} = d_{4,0} \ge 2r_g$ and Appendix D). In the case where $d_{2,0} = d_{4,0} \ge 2r_g$, the connectivity of forces within the cell is maintained through contact with neighbouring cells. This is reflected by the external force $F_1^a$ being balanced by the tangential external force $G_2$ (see Fig. 2), resulting in the equilibrium of the two half cells (see Appendix D). As a result, the initial density in the experiments can be accounted for by calibrating the two numerical parameters $d_{2,0}$ and $\alpha_0$. The sintering time before loading is included in the snow-H-model from the initial bond radius $r_{b,0}$, which is capped by the grain radius. The numerical parameters and their range of variation are summarized in the last two columns of Table 1.

The first objective of the calibration step is to determine the initial value of the three numerical parameters describing both the geometry and the sintering state of the H-cell, namely the opening angle $\alpha_0$, the intergranular distance $d_{2,0}$ and the bond radius $r_{b,0}$. To this end, stress-strain curves for different initial parameters are systematically compared to the experimental curves by Abele and Gow (1976). The best fitting triplet (denoted by the index bf) in the parameter set $\Omega$ given in Table 1 corresponds to the triplet $(\alpha_0, d_{2,0}, r_{b,0})_{bf}$ which verifies:

$$ \text{err}(\alpha_0, d_{2,0}, r_{b,0})_{\text{bf}} = \min_{\Omega} (\text{err}(\alpha_0, d_{2,0}, r_{b,0})) \quad (32) $$

with err being a relative error gap function defined by:

$$ \text{err} = \frac{\sum_{i=1}^{n_{\text{exp}}} |\sigma_{\text{num}}(\varepsilon_i) - \sigma_{\text{exp}}(\varepsilon_i)| (\varepsilon_i - \varepsilon_{i-1})}{\sum_{i=1}^{n_{\text{exp}}} |\sigma_{\text{exp}}(\varepsilon_i)| (\varepsilon_i - \varepsilon_{i-1})} \quad (33) $$

In Eqn (33), the numerator corresponds to the cumulated area between the numerical and experimental curves in the axial stress-strain space $(\sigma, \varepsilon)$, whereas the denominator corresponds to the area below the experimental curve. $n_{\text{exp}}$ is the number of points taken from the experimental stress-strain curve, $\varepsilon_i$ is the strain at the $i$th point of the experimental curve, $\sigma_{\text{num}}$ is the stress obtained with the 3-D H-model and $\sigma_{\text{exp}}$ is the experimental one.

The parameter set in $\Omega = \{(\alpha_0, d_{2,0}, r_{b0}) \in [45^{\circ}; 90^{\circ}] \times [2r_g; 3.6r_g] \times [0.05r_g; 0.9r_g]\}$ has been explored for experimental curves corresponding to a density between 470 kg/m$^3$ and 640 kg/m$^3$. The stress-strain curves of the experimental data and the corresponding best fitting numerical curves in $\Omega$ are reported in Fig. 7.

Table 2 reports the best fit parameters with the initial density of the experimental specimens increasing from 470 kg/m$^3$ to

640 kg/m$^3$. The relative error gap, as defined in Eqn (33), ranges between 1.4% and 12.3%.

When plotted as a function of the initial snow density in Fig. 7, the error is found smaller in the middle of the density range. In any case, the value of the relative error between experimental and numerical curves remains relatively small, which confirms the ability of the 3-D H-model to be used as a constitutive model for snow, at least within the density range [470; 640] kg/m$^3$. Higher densities are impossible to reach with the 3-D H-model, as the condition $\alpha_0 = 45^{\circ}$ and $d_{2,0}/r_g = 2$ corresponds to the limit case where a contact between the two hexagons within a cell is created. For lower densities, the mismatch between experimental and numerical curves becomes too large to include these densities in the range of validity of the snow-H-model.

This first calibration was carried out blindly, in the sense that only the minimization of the gap between measured and predicted data was targeted, without considering the agreement between experimental density and model density, the latter being controlled by the two geometric parameters of the H-cell. This is the purpose of the next section, to achieve a much more versatile model.

## 4. Discussion on the model versatility

As the snow-H-model has demonstrated its ability to reproduce experimental results through parameter adjustment, the next step is to assess its predictive capability for snow behaviour under varying conditions of density, temperature and strain rate.

### 4.1. Relation between mesoscopic geometry and macroscopic density

Calculating the density in multi-scale models is challenging as its mesostructure only includes those elements which are actively participating in the force transmission. For instance, in the present model, the bi-hexagonal pattern is designed as the minimal structure to depict only snow grains which are part of force chains. The surrounding space includes both void and snow grains not participating to force transmission. This is the reason why we cannot and should not calibrate the H-cell geometry to match directly the macroscopic density. Therefore, we have chosen first to calibrate the initial geometrical parameters of the cells to match the experimental stress-strain curves from Abele and Gow (1976). The objective of this section is to provide an estimation of the proportion of ice grains participating to force transmission by comparing the macroscopic snow density $\rho_{\text{snow}}$ and the mesoscopic density $\rho_{\text{meso}}$.

In the snow-H-model, we selected the smallest peripheral parallelepiped volume as the H-cell volume $V_{\text{meso}}$ (see Fig. 9). Under this assumption, the mesoscopic density can be expressed as follows:

$$ \rho_{\text{meso}}(\alpha_0, d_{2,0}) = \frac{V_{\text{ice}}}{V_{\text{meso}}(\alpha, d_{2,0})} \rho_{\text{ice}} = \frac{40\pi r_g^3}{3L_1L_2L_3} \rho_{\text{ice}} \quad (34) $$

where $V_{\text{ice}}$ is the volume of solid ice, which corresponds to the volume of the ten spherical grains (we neglect the ice bond volumes), and $\rho_{\text{ice}}$ is the density of ice, equal to 920 kg/m$^3$.

The best fitting parameters obtained in Section 3.2 provide an initial density $\rho_{\text{meso}}$ lower than 250 kg/m$^3$, while the initial densities measured in the experiments are larger than 470 kg/m$^3$.

The density based on the bounding box shown in Fig. 8 systematically underestimates the real snow density for several reasons:
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Stress-strain curves for confined compression test: experimental results from Abele and Gow (1976) and best fitting numerical curves obtained with the 3-D H-model.</figcaption></figure>

<table><caption>**Table 2.** Parameters giving the best fit for each experimental curve in $\Omega = [45^{\circ}; 90^{\circ}] \times [2 r_g; 3.6 r_g] \times [0.05 r_g; 0.9 r_g]$.</caption><thead><tr><th rowspan="2">Exp Nb*</th><th colspan="4">Data from Abele and Gow (1976)</th><th colspan="4">Best fit parameters for S ∈ Ω</th></tr><tr><th>Density (kg m<sup>-3</sup>)*</th><th>Storage temperature (°C)*</th><th>Strain rate (s<sup>-1</sup>)*</th><th>Snow age (days)*</th><th>Initial opening angle (°)</th><th>Intergranular distance (/r<sub>g</sub>)</th><th>Initial bond radius [0; r<sub>g</sub>]</th><th>Gap error (%)</th></tr></thead><tbody><tr><td>50</td><td>470</td><td>-1</td><td>7.84</td><td>3</td><td>51.0</td><td>2.75</td><td>0.50</td><td>3.3</td></tr><tr><td>48</td><td>480</td><td>-34</td><td>7.84</td><td>3</td><td>53.0</td><td>3.05</td><td>0.03</td><td>10.5</td></tr><tr><td>52</td><td>490</td><td>-1</td><td>7.84</td><td>3</td><td>63.5</td><td>3.15</td><td>0.40</td><td>2.4</td></tr><tr><td>51</td><td>520</td><td>-1</td><td>-7.84</td><td>3</td><td>50.0</td><td>2.90</td><td>0.45</td><td>1.4</td></tr><tr><td>35</td><td>530</td><td>-34</td><td>5.26</td><td>7</td><td>51.5</td><td>2.85</td><td>0.30</td><td>4.2</td></tr><tr><td>41</td><td>580</td><td>-1</td><td>16.0</td><td>7</td><td>55.0</td><td>2.10</td><td>0.65</td><td>2.1</td></tr><tr><td>36</td><td>640</td><td>-34</td><td>5.26</td><td>7</td><td>45.0</td><td>2.00</td><td>0.60</td><td>12.3</td></tr></tbody></table>

<figure><img src="image_10.png" /><figcaption>**Figure 8.** Gap error between experimental curves and best fitting numerical curves with initial parameters not respecting the density condition (blue crosses) and respecting the density condition presented in Section 4.1 (orange circles) as a function of the initial experimental density.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 9.** Definition of the peripheral volume of the H-cell $V_{meso}$.</figcaption></figure>

1. As already mentioned, the H-cell accounts only for grains participating in stress transmission (the force chains) whereas, in real granular materials, many grains do not transmit significant stresses. Therefore, many ice grains are not accounted for by the bi-hexagonal mesostructure, thus underestimating the solid volume.

2. The arbitrary choice for the mesoscopic volume includes much empty space and probably overestimates the void volume.

3. To a lower extent, the ideal spherical shape of the ice grain in the H-model approximates the shape of the real ice grains. In real snow, ice grains are neither spherical and nor

monodisperse (see, for instance, Gay and others (2002), Kaempfer and Schneebeli (2007) and Calonne and others (2014)).

To account for the additional ice in the pore space, a coefficient $\gamma$ can be introduced in the relationship between macroscopic and mesoscopic densities, as:

$$ \rho_{\text{macro}} = \gamma \rho_{\text{meso}}. \qquad (35) $$

<footnote>*Data from Abele and Gow (1976) experiments.</footnote>
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Density as a function of the opening angle for different initial intergranular distance $d_2$. Lines represent the analytical expression of the macroscopic density with $\gamma = 2.58$. Each symbol represents a simulation, where the y-coordinate is the snow density of the corresponding experimental tests.</figcaption></figure>

<table><caption>**Table 3.** Best fitting parameters over $S_{bf,c}$ with the additional density constraint for the same experiment from Abele and Gow (1976) used in Section 3</caption><thead><tr><th colspan="5">Data from Abele and Gow (1976)</th><th colspan="4">Best fit parameters with density constraint</th></tr><tr><th>Exp Nb</th><th>Density (kg m<sup>-3</sup>)</th><th>Temperature (°C)</th><th>Strain rate (s<sup>-1</sup>)</th><th>Snow age (days)</th><th>Initial opening angle (°)</th><th>Intergranular distance (/r<sub>g</sub>)</th><th>Initial bond radius (/r<sub>g</sub>)</th><th>Gap error (%)</th></tr></thead><tbody><tr><td>50</td><td>470</td><td>-1</td><td>7.84</td><td>3</td><td>73.5</td><td>3.1</td><td>0.01</td><td>13.9</td></tr><tr><td>48</td><td>480</td><td>-34</td><td>7.84</td><td>3</td><td>53.0</td><td>3.3</td><td>0.35</td><td>12.6</td></tr><tr><td>52</td><td>490</td><td>-1</td><td>7.84</td><td>3</td><td>48.0</td><td>3.5</td><td>0.65</td><td>2.7</td></tr><tr><td>51</td><td>520</td><td>-1</td><td>7.84</td><td>3</td><td>50.0</td><td>2.9</td><td>0.45</td><td>1.5</td></tr><tr><td>35</td><td>530</td><td>-34</td><td>5.26</td><td>7</td><td>51.5</td><td>2.6</td><td>0.20</td><td>10.8</td></tr><tr><td>41</td><td>580</td><td>-1</td><td>16.0</td><td>7</td><td>45.0</td><td>2.5</td><td>0.35</td><td>7.3</td></tr><tr><td>36</td><td>640</td><td>-34</td><td>5.26</td><td>7</td><td>45.0</td><td>2.0</td><td>0.6</td><td>12.3</td></tr></tbody></table>

This coefficient is set to $\gamma = 2.58$ from the harmonic average of the ratios $\frac{\rho_{\text{exp}}}{\rho_{\text{meso}}}$ based on the calibration tests shown in Section 3 (see Table 2):

$$ \gamma = \left( \frac{1}{n_{\text{test}}} \sum_{i=1}^{n_{\text{test}}} \frac{\rho_{\text{meso},i}}{\rho_{\text{exp},i}} \right)^{-1} \quad (36) $$

Figure 10 displays the analytical expression of the macroscopic density with $\gamma = 2.58$, as a function of the opening angle for different intergranular distances, with each test used in previous Section 3.2 corresponding to a marker. The relative error is less than 2% for each test.

The high value of $\gamma$ (with respect to 1) confirms that a considerable percentage of ice material is not accounted for in the density calculation. Specifically, $\gamma = 2.58$ implies that the H-cells account for 38% of the mass of the material. This order of magnitude almost aligns with the observations by Wautier and others (2017) that the fraction of grains contributing to force chains ranges between 23% and 31% in sand.

In contrast to Section 3.2, the fitting procedure can be optimized by adding the additional constraint that the model density matches the experimental one. Thus, for each test i, the best fitting triplet of initial geometrical parameters that satisfies $\rho_{\text{macro}}(\alpha_0, d_{2,0}) = \rho_{\text{exp},i}$ can be searched in $\Omega_c(\rho)$ defined by:

$$ \Omega_c(\rho) = \{(\alpha_0, d_{2,0}, r_{b,0}) \in \Omega | \rho_{\text{macro}}(\alpha_0, d_{2,0}) = \rho\} \quad (37) $$

The best fitting parameters accounting for the density relationship for each test are summarized in Table 3 and the corresponding

stress-strain curves are plotted in Fig. 11. The maximum relative error observed is of 13.9%, which is only slightly larger than without the constraint (39) (see Fig. 7). In general, the largest errors correspond to the extremal values of densities, while low gap errors are found for $\rho_{\text{exp}}$ close to 500 kg/m$^3$. The discrepancy between experimental and numerical curves is maximal for small and large deformations. At small deformation, an obvious inflexion point is observed in numerical curves at the transition between elastic and plastic regimes, due to the breakage of ice bonds. Indeed, as all bridges share the same radius, ice bridges all break at the same time for a given direction of mesostructure. In a real snow sample, the radius of ice bridges varies, which means they do not break at the same deformation level. Consequently, a smoother transition between elastic and plastic regime is observed experimentally. For large deformation, the H-model does not account for the reorganization of the mesostructure, as the cells are deformed independently. This limitation of the model has been discussed by Deng (2022), who introduced a way to reinitialize mesostructure geometries to increase the domain of application.

## 4.2. Relative influence of the three microstructure parameters on the mechanical response

Among the three parameters of the snow-H-model, two of them control the initial density (namely the initial intergranular distance $d_{2,0}$ and the initial opening angle $\alpha_0$). The remaining parameter is
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Stress-strain curves for confined compression tests: experimental results from Abele and Gow (1976) and best fitting numerical curves obtained with the 3-D H-model for initial parameters verifying $\rho_{\text{macro}}(\alpha_0, d_{2,0}) = \rho_{\text{exp}}$.</figcaption></figure>

the initial bond radius $r_{b,0}$. It accounts for the impact of temperature and age of the snow sample. In this subsection, we characterize the relative influence of the three initial parameters of the 3-D H-model on the stress-strain curve in oedometer conditions.

## 1. Impact of the initial density

The initial density of a numerical snow sample is primarily determined by the geometry of the H-cells, characterized by the parameters $\alpha_0$ and $d_{2,0}$. There is a direct relationship between $d_2$ and density, as the void within an H-cell increases proportionally with $d_2$. The stress-strain curves for a confined compression test with different initial intergranular distances $d_{2,0}$ are plotted in Fig. 12 (up) while the two other parameters are kept constant. Increasing the intergranular distance reduces the mesoscopic density and results in smaller axial stresses (as shown in Fig. 12 (up)), as it was observed when reducing the density of a snow sample (Landauer, 1955; Radke and Hobbs, 1967; Abele and Gow, 1975; Scapozza and Bartelt, 2003; Wautier and others, 2015).

On the other hand, the relationship between the opening angle $\alpha$ and the density is not bijective. For small values of $\alpha$, density increases with decreasing $\alpha$; for larger values of $\alpha$, density increases with $\alpha$ (see Fig. 10). In Fig. 12 (down), a compression test with $d_2/r_g = 3.5$ and $r_{b,0}/r_g = 0.65$ was run with different initial opening angles $\alpha_0$. The evolution of axial stress with opening angle follows the same pattern as for density. For opening angles lower than 55°, axial stress decreases when the opening angle is increased while the opposite trend is observed for opening angles larger than 55°. However, the mesoscopic density increases with the opening angle for opening angles significantly larger than 55° (between 60° and 70°, depending on the initial intergranular distance $d_{2,0}$). Thus, the non-monotonic influence of the opening angle on the stress-strain curves probably reflects the complex influence of density and other microstructural parameters (e.g. anisotropy of the contact distribution).

## 2. Impact of the sintering time

Another key geometrical parameter to calibrate is the initial bond radius $r_{b,0}$. While the sintering process that occurs prior to mechanical loading is not directly modelled, the initial bond radius serves as a proxy for this process since it has been shown that bond

radius grows over sintering time (Kuroiwa, 1961; Herwijnen and Miller, 2013), leading to stronger bonds (Szabo and Schneebeli, 2007). In Fig. 13, the same compression test is now considered with different initial bond radii. It can be observed that, at the beginning of the loading, the initial bond radius has a positive effect on the axial stress. Then, as ice bonds are broken under compression, this influence of the initial bond radius is progressively reduced. It should be emphasized that the observation of enhanced stresses with the initial bond radius is somehow reminiscent to snow stiffness evolution during isothermal metamorphism (Wautier and others, 2015).

Note that the growth in bond radius with sintering time cannot be readily observed from the data reported in Table 3, where the values of $r_{b,0}$ for 3 day and 7 day samples are not significantly different. This can be attributed to the fact that in the experimental work of Abele and Gow (1976), the snow samples were not collected immediately after a snowfall and sintering occurred prior to the sample collection. Consequently, it is difficult to link the degree of sintering solely to the storage duration.

## 4.3. Influence of loading conditions on the mechanical response

In this subsection, we analyse the ability of the snow-H-model to account for loading conditions that differ from those of Abele and Gow (1976). We focus on the strain rate and temperature influence on the mechanical response in oedometer conditions.

### 1. Impact of the strain rate

As the ice bonds are described by a visco-elastoplastic model, it is relevant to investigate the effect of the loading rate on the material response. As already specified, the analysis considers exclusively loading rates above $10^{-4} \text{s}^{-1}$ to remain within the brittle regime (Schulson, 2001). We can define the characteristic time $t_{\text{viscous}}$ at which the orders of magnitude of the viscoplastic and elastic deformations coincide:

$$t_{\text{viscous}} = \frac{\varepsilon_{b,n}^{e}}{A(T_{\text{emp}}) \sigma_{b,n}^{m}} \approx \frac{\sigma_{\text{bcs}}}{E_{\text{oedo}} A(T_{\text{emp}}) \sigma_{\text{bcs}}^{m}} \quad (38)$$
<figure><img src="image_14.png" /><figcaption>**Figure 12.** Stress-strain curves for confined compression test with different initial intergranular distances, $\alpha_0 = 48^\circ$ and $r_{b,0}/r_g = 0.65$ (up), and with different initial opening angles, $d_{2,0}/r_g = 3.5$ and $r_{b,0}/r_g = 0.65$ (down).</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 13.** Stress-strain curves for confined compression test with different initial bond radius, $\alpha_0 = 48^\circ$ and $d_{2,0}/r_g = 3.5$.</figcaption></figure>

where $\sigma_{bcs}$ is the brittle compression strength of ice. At $-10^\circ\text{C}$, Eqn (38) gives $t_{\text{viscous}} = 24$ s with $E_{\text{oedo}} = 1.6$ GPa the calibrated value of oedometer coefficient, and $\sigma_{\text{bcs}} = 13.75$ MPa from Eqn (27).

In the experiment by Abele and Gow (1976), the strain rates ranged between $5\text{s}^{-1}$ and $16\text{s}^{-1}$, resulting in a total loading time less than 0.08 s, thus much smaller than $t_{\text{viscous}}$. Consequently, for such high strain rates, the effects of viscosity can be considered negligible.

The stress-strain curves of confined compression tests performed with different strain rates on similar numerical snow samples ($\alpha_0 = 48^\circ$, $r_{b,0} = 0.6r_g$, $d_{2,0} = 3.54r_g$) are given in Fig. 14. To measure the direct influence of viscosity, these curves are compared with a snow model with nonviscous bonds. The evolution of the percentage of broken bonds during the loading for the different strain rates is reported in Fig. 14. The axial stress increases

with the strain rate up to $\dot{\varepsilon} = 0.025\text{s}^{-1}$. For higher values, the strain rate has no significant influence on the stress-strain response which is the same as the response of the nonviscous model. This is consistent with the fact that the characteristic time of the viscous deformation $t_{\text{viscous}}$ is larger than the time needed to break the bonds at high strain rates: for example, with $\dot{\varepsilon} = 0.01\text{s}^{-1}$, more than 70% of the bonds are broken after only 3 s, as shown in Fig. 14. The proportion of broken bonds increases with the strain rate (Schulson, 2001), as the loading rate controls the stress evolution within the bond through Glen's law (Eqn (19)). As the contact between unbonded grains is assumed to be rate-independent in the proposed model, the strain rate effect decreases as the number of bond breakage rises, bringing the curves to converge at high strain rates. It is worth noting that at the lowest strain rates considered, the test duration could allow for the formation of additional bonds at the new contact points (see Appendix C). While the current model does not account for the formation of these new bonds, their inclusion would probably result in a more significant impact of viscosity at low strain rates.

## 2. Impact of the temperature

The second loading parameter to studied is the temperature. McClung (1996), Schweizer (1998) and Takei and Maeno (2004) all observed that an increase in temperature leads to a decrease in snow strength and stiffness. The temperature affects the snow mechanical behaviour independently during the sintering phase and the loading phase.

During the sintering phase, it was shown that a higher temperature leads to a higher bond growth rate (Blackford, 2007; Herwijnen and Miller, 2013) resulting in a greater sintering force
<figure><img src="image_16.png" /><figcaption>**Figure 14.** Stress-strain curves (left) and evolution of the proportion of broken bonds (right) for confined compression tests at different strain rates, compared with a snow material with nonviscous bonds.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Figure 15.** Stress-strain curves (left) and evolution of the proportion of broken bonds (right) for confined compression tests at different temperatures and for $\dot{\epsilon} = 7.84 \text{ s}^{-1}$.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Figure 16.** Stress-strain curves for confined compression tests at different temperatures and strain rates.</figcaption></figure>

and to an increase in the bond radius (Szabo and Schneebeli, 2007; Bahaloo, 2022). As the snow-H-model does not model the sintering phase, the effect of the temperature during sintering cannot be assessed. However, the initial bond radius can be used as a proxy. We deduce from Fig. 13 that we have an increase of the stress with the temperature during the sintering time.

During the mechanical loading, the temperature influences both the brittle compression strength of ice defined by Eqn (28) (Schulson, 2001) and the viscous behaviour of the ice bonds through Glen's law (Eqn (19)).

To measure the influence of the temperature in the model, isotropic compression tests were simulated from the same numerical snow sample ($\alpha_0 = 48^\circ$, $r_{b,0} = 0.6r_g$, $d_{2,0} = 3.54r_g$) with a strain rate of $\dot{\epsilon} = 7.84 \text{s}^{-1}$ at different temperatures (Fig. 15). It can be consistently observed that, for a given axial strain, the axial stresses are larger for lower temperatures, as bond failure occurs earlier in warmer conditions, as illustrated in Fig. 15 and described in Schulson (1990).

Finally, Fig. 16 shows the stress-strain curves for different temperatures with a strain rate of $7.84 \text{s}^{-1}$, $0.001 \text{s}^{-1}$ and $0.0001 \text{s}^{-1}$.
The strain rate has no significant effect for temperatures strictly below −10°C. At higher temperatures, namely $T_{\text{emp}} = -1$°C and $T_{\text{emp}} = -10$°C, the effect of the strain rate becomes apparent. This stems from the temperature dependent coefficient $A(T_{\text{emp}})$ in Glen's law (Eqn (19)). As this coefficient increases with temperature, the characteristic time of the viscous deformation evolves inversely with temperature. As previously mentioned, the discrepancy between the curves at $T_{\text{emp}} = -1$°C and $T_{\text{emp}} = -10$°C is larger at the lowest strain rate, as the effect of loading rate on ice behaviour is more marked.

From Figs. 15 and 16, we can conclude that, for a given axial strain, the axial stress increases with decreasing temperature, whereas the snow stiffness is not affected by temperature. However, the expected influence of temperature on snow stiffness could be reproduced by considering in the model the known temperature-dependency of ice stiffness (Parameswaran, 1987).

## 5. Conclusions and perspectives

The 3-D H-model has been extended to snow by adding an ice bond component to the usual frictional contact law. A series of confined compression tests have been performed from numerical snow samples. The axial stress results have been compared against the experimental curves of Abele and Gow (1976). The ability of the snow-H-model to reproduce the mechanical behaviour of snow has been demonstrated for the intermediate range of snow density [470; 640] kg/m³, at different temperature and strain rates. For snow of lower or larger density, the calibration of the model is no longer adequate as, constitutively, the snow H-model is only able to represent the snow microstructure of intermediate density. To describe both loose and dense snow specimens, the H-cell geometry would need further adaptations. For example, changing the number of grains in the H-cell could allow to reach higher or lower density. It is worth noting that within the considered density range, the assumption of no rotation of the ice grains is reasonable, as increasing the material density enhances the connectivity between grains.

Within the domain of validity of the snow-H-model, the relationship between the geometrical parameters of the H-cell and the macroscopic density has been studied closely. Indeed, the statistical nature of the 3-D H-model prevents the pore space between the different H-cells from being accounted for. Consequently, the relationship between the mesoscopic and the macroscopic densities was derived using a proportionality coefficient. This coefficient accounts for the ice fraction that does not contribute to the load-bearing capacity of snow (i.e. the ice grains not participating to the force chains) and the void in between the H-cells.

After calibration of the model against the experimental work of Abele and Gow (1976), the effect of the initial bond radius has been investigated, revealing that a change in the initial bond radius produces a phenomenological response similar to what would be observed with a change in the aging time of a sample under isothermal conditions.

Finally, the impact of the loading conditions (strain rate and temperature) on the response of the snow-H-model has been highlighted. It has been shown that the strain rate (when sufficiently low) directs a substantial effect on both the viscous response of the bonds and the bond failure rate. Increasing temperature also leads to faster bond failure and smaller stresses in the material. It should be noted that the influence of strain rate on the stress-strain response is greater at higher temperatures.

As the capacity of the snow-H-model to reproduce experimental oedometer test has been proven, the next step will consist in studying the predictive ability of the model. This would require additional experiments using snow sample with a better control of the initial microstructure, and a larger range of loading condition, with lower strain rates. An implementation of the extended version of the 3-D H-model in a computational software suitable for dealing with practical engineering concerns is also considered. Using such a multi-scale approach ranging from the particle to the continuum scale is thought to bridge the gap between engineering-scale systems and the microscopic features of snow. Comparing the H-model with DEM simulations at a smaller intermediate scale would help assess computational efficiency and accuracy for snow model, as has been already demonstrated for sand on large geotechnical systems involving scales beyond the reach of DEM simulations (Xiong and others, 2019).

The multi-scale approach of the model can allow to study the effect of snow microstructure in more detail. Here, the effects of the microparameters related to the snow density have been widely studied, but it would be also possible to study the effect of anisotropy. All the simulations in this paper have been realized with an isotropic distribution of direction, but it would be possible to include anisotropy by modifying the probability distribution function $\omega(\theta, \varphi, \psi)$.

New extensions of the snow-H-model can also be considered in the future. Ice bond creation and ductile failure could be added in the model, to extend the range of applicability of the model to lower strain rate. Such an improvement would make it possible to consider, for example, avalanche problems and to cover more domains application of snow mechanics.

**Acknowledgements.** This work was supported by the Luxembourg National Research Foundation FNR (SNOWTEC, project ref. 16241613).

A CC-BY public copyright license has been applied by the authors to the present document and will be applied to all subsequent versions up to the Author Accepted Manuscript arising from this submission, in accordance with the grant's open access conditions.

The authors would like to thank the anonymous reviewers for their deep reviews of the manuscript. We sincerely appreciate the intensive work done, which helped us improve the quality of the manuscript.

**Competing interests.** The authors declare no competing interests.

## References

* Abele G and Gow A (1975) Compressibility characteristics of undisturbed snow. https://apps.dtic.mil/sti/tr/pdf/ADA012113.pdf (accessed 31 January 2025).

* Abele G and Gow AJ (1976) *Compressibility Characteristics of Compacted Snow*. Hanover, N.H: CRREL.

* Bahaloo H and 6 others (2022) Ice sintering: Dependence of sintering force on temperature, load, duration, and particle size. *Journal of Applied Physics* 131(2), 025109. doi: 10.1063/5.0073824

* Barnes P, Tabor D and Walker JCF (1971) The friction and creep of polycrystalline ice. *Proceedings of the Royal Society of London A Mathematical and Physical Sciences* 324(1557), 127–155. doi: 10.1098/rspa.1971.0132

* Blackford JR (2007) Sintering and microstructure of ice: A review. *Journal of Physics D Applied Physics* 40(21), R355–R385. doi: 10.1088/0022-3727/40/21/r02

* Blatny L, Löwe H and Gaume J (2023) Microstructural controls on the plastic consolidation of porous brittle solids. *Acta Materialia* 250, 118861. doi: 10.1016/j.actamat.2023.118861

* Blatny L, Löwe H, Wang S and Gaume J (2021) Computational micromechanics of porous brittle solids. *Computers and Geotechnics* 140, 104284. doi: 10.1016/j.compgeo.2021.104284
Bobillier G, Bergfeld B, Dual J, Gaume J, van Herwijnen A and Schweizer J (2021) Micro-mechanical insights into the dynamics of crack propagation in snow fracture experiments. *Scientific Reports* 11(1), 11711. doi: 10.1038/s41598-021-90910-3

Calonne N, Flin F, Geindreau C, Lesaffre B and du Roscoat SR (2014) Study of a temperature gradient metamorphism of snow from 3-D images: Time evolution of microstructures, physical properties and their associated anisotropy. *The Cryosphere* 8, 2255–2274.

Cundall PA and Strack ODL (1979) A discrete numerical model for granular assemblies. *Géotechnique* 29(1), 47–65. doi: 10.1680/geot.1979.29.1.47

Deng N (2022) Micromechanics of critical state and its emergence in multiscale modeling of soil failure. PhD Thesis, Universite Grenoble Alpes.

Dent JD and Lang TE (1980) Modeling of snow flow. *Journal of Glaciology* 26(94), 131–140. doi: 10.3189/S0022143000010674

El-Sayegh Z and El-Gindy M (2019) Modelling and prediction of tyre–snow interaction using finite element analysis–smoothed particle hydrodynamics techniques. *Proceedings of the Institution of Mechanical Engineers, Part D: Journal of Automobile Engineering* 233(7), 1783–1792. doi: 10.1177/0954407018788997

Fauve M, Rhyner H, Lüthi A, Schneebeli M and Lehning M (2008) Putting snow knowledge into the development of winter sports equipment. *Sports Technology* 1(2–3), 145–151. doi: 10.1080/19346182.2008.9648465

Fierz C and 8 others (2009) The international classification for seasonal snow on the ground. UNESCO, IHP (International Hydrological Programme)–VII, Technical Documents in Hydrology, No 83; IACS (International Association of Cryospheric Sciences) contribution No 1.

Fourteau K, Freitag J, Malinen M and Löwe H (2024) Microstructure-based simulations of the viscous densification of snow and firn. *The Cryosphere* 18(6), 2831–2846. doi: 10.5194/tc-18-2831-2024

Gaume J, Chambon G, Eckert N, Naaim M and Schweizer J (2015) Influence of weak layer heterogeneity and slab properties on slab tensile failure propensity and avalanche release area. *The Cryosphere* 9(2), 795–804. doi: 10.5194/tc-9-795-2015

Gaume J, Chambon G, Naaim M and Eckert N (2011) Influence of weak layer heterogeneity on slab avalanche release using a finite element method. In Bonelli S, Dascalu C and Nicot F (eds), *Advances in Bifurcation and Degradation in Geomaterials*. Netherlands, Dordrecht: Springer, 261–266.

Gaume J, Gast T, Teran J, van Herwijnen A and Jiang C (2018) Dynamic anti-crack propagation in snow. *Nature Communications* 9(1), 3047. doi: 10.1038/s41467-018-05181-w

Gaume J, van Herwijnen A, Chambon G, Birkeland KW and Schweizer J (2015) Modeling of crack propagation in weak snowpack layers using the discrete element method. *The Cryosphere* 9(5), 1915–1932. doi: 10.5194/tc-9-1915-2015

Gaume J, van Herwijnen A, Gast T, Teran J and Jiang C (2019) Investigating the release and flow of snow avalanches at the slope-scale using a unified model based on the material point method. *Cold Regions Science and Technology* 168, 102847. doi: 10.1016/j.coldregions.2019.102847

Gay M, Fily M, Genthon C, Frezzotti M, Oerter H and Winther J-G (2002) Snow grain-size measurements in Antarctica. *Journal of Glaciology* 48, 527–535.

Guillet L, Blatny L, Trottet B, Steffen D and Gaume J (2023) A depth-averaged material point method for shallow landslides: Applications to snow slab avalanche release. *Journal of Geophysical Research: Earth Surface* 128(8), e2023JF007092. doi: 10.1029/2023JF007092

Haehnel RB and Shoop SA (2004) A macroscale model for low density snow subjected to rapid loading. *Cold Regions Science and Technology* 40(3), 193–211. doi: 10.1016/j.coldregions.2004.08.001

Herwijnen AV and Miller DA (2013) Experimental and numerical investigation of the sintering rate of snow. *Journal of Glaciology* 59(214), 269–274. doi: 10.3189/2013jog12j094

Hjort M, Eriksson O and Bruzelius F (2017) Comprehensive study of the performance of winter tires on ice, snow, and asphalt roads: The influence of tire type and wear. *Tire Science and Technology* 45, 175–199. doi: 10.2346/tire.17.450304

Johnson JB and Hopkins MA (2005) Identifying microstructural deformation mechanisms in snow using discrete-element modeling. *Journal of Glaciology* 51(174), 432–442. doi: 10.3189/17275650578129188

Kabore BW, Peters B, Michael M and Nicot F (2021) A discrete element framework for modeling the mechanical behaviour of snow—Part I: Mechanical behaviour and numerical model. *Granular Matter* 23(2), 42. doi: 10.1007/s10035-020-01083-1

Kaempfer TU and Schneebeli M (2007) Observation of isothermal metamorphism of new snow and interpretation as a sintering process. *Journal of Geophysical Research: Atmospheres* 112(D24). doi: 10.1029/2007JD009047

Kojima K (1967) Densification of seasonal snow cover. *Physics of Snow and Ice: Proceedings* 1(2), 929–952.

Kruyt NP and Rothenburg L (1996) Micromechanical definition of the strain tensor for granular materials. *Journal of Applied Mechanics* 63(3), 706–711. doi: 10.1115/1.2823353

Kuroiwa D (1961) A study of ice sintering. *Tellus A: Dynamic Meteorology and Oceanography* 13(2), 252–259. doi: 10.3402/tellusa.v13i2.9450

Landauer JK (1955) Stress-strain relations in snow under uniaxial compression. *Journal of Applied Physics* 26(12), 1493–1497. doi: 10.1063/1.1721937

Lee JH (2009) A new indentation model for snow. *Journal of Terramechanics* 46(1), 1–13. doi: 10.1016/j.jterra.2009.02.001

Lee JH (2013) Calibration and validation of a tire–snow interaction model. *Journal of Terramechanics* 50(5), 289–302. doi: 10.1016/j.jterra.2013.09.003

Lipovsky BP (2022) Density matters: Ice compressibility and glacier mass estimation. *Journal of Glaciology* 68(270), 831–832. doi: 10.1017/jog.2021.132

Liu J, Nicot F and Zhou W (2018) Sustainability of internal structures during shear band forming in 2D granular materials. *Powder Technology* 338, 458–470.

Liu Z, Nicot F, Wautier A and Darve F (2024) Multiscale investigation of bonded granular materials: The H-bond model. *Computers and Geotechnics* 172, 106481. doi: 10.1016/j.compgeo.2024.106481

Love A (1892) A treatise on the mathematical theory of elasticity. https://hal.archives-ouvertes.fr/hal-01307751 (accessed 31 January 2025).

Luding S (2008) Introduction to discrete element methods. *European Journal of Environmental and Civil Engineering* 12(7–8), 785–826. doi: 10.1080/19648189.2008.9693050

Ma Q, Wautier A and Nicot F (2022) Mesoscale investigation of fine grain contribution to contact stress in granular materials. *Journal of Engineering Mechanics* 148(3), 04022005. doi: 10.1061/(ASCE)EM.1943-7889.0002078

McClung DM (1996) Effects of temperature on fracture in dry slab avalanche release. *Journal of Geophysical Research: Solid Earth* 101(B10), 21907–21920. doi: 10.1029/95JB03114

Mede T, Chambon G, Hagenmuller P and Nicot F (2018) Snow failure modes under mixed loading. *Geophysical Research Letters* 45(24), 13,351–13,358. doi: 10.1029/2018GL080637

Meschke G, Liu C and Mang HA (1996) Large strain finite-element analysis of snow. *Journal of Engineering Mechanics* 122(7), 591–602. doi: 10.1061/(ASCE)0733-9399(1996)122:7(591)

Michael M (2014) A discrete approach to describe the kinematics between snow and a tire tread. PhD Thesis, Universite du Luxembourg, The Faculty of Sciences, Technology and Communication.

Mulak D and Gaume J (2019) Numerical investigation of the mixed-mode failure of snow. *Computational Particle Mechanics* 6(3), 439–447. doi: 10.1007/s40571-019-00224-5

Nicot F (2004) From constitutive modelling of a snow cover to the design of flexible protective structures Part I—Mechanical modelling. *International Journal of Solids and Structures* 41(11), 3317–3337. doi: 10.1016/j.ijsolstr.2003.09.022

Nicot F and Darve F (2011) The H-microdirectional model: Accounting for a mesoscopic scale. *Mechanics of Materials* 43(12), 918–929. doi: 10.1016/j.mechmat.2011.07.006

Parameswaran VR (1987) Orientation dependence of elastic constants for ice. *Defence Science Journal* 37(3), 367–375. doi: 10.14429/dsj.37.5924

Peters B, Kabore BW, Michael M and Nicot F (2021) A discrete element framework for modeling the mechanical behaviour of snow Part II: Model validation. *Granular Matter* 23(2), 43. doi: 10.1007/s10035-020-01084-0

Podolskiy E, Chambon G, Naaim M and Gaume J (2013) A review of finite-element modelling in snow mechanics. *Journal of Glaciology* 59, 1189–1201. doi: 10.3189/2013JoG13J121
Radjai F, Wolf DE, Jean M and Moreau J-J (1998) Bimodal character of stress transmission in granular packings. *Physical Review Letters* 80(1), 61–64. doi: 10.1103/PhysRevLett.80.61

Radke LF and Hobbs PV (1967) The strength–density relationship for dry snow. *Journal of Glaciology* 6(48), 893–896. doi: 10.3189/S0022143000020 189

Satake M (1992) A discrete-mechanical approach to granular materials. *International Journal of Engineering Science* 30(10), 1525–1533. doi: 10.1016/ 0020-7225(92)90162-A

Scapozza C and Bartelt P (2003) Triaxial tests on snow at low strain rate. Part II. Constitutive behaviour. *Journal of Glaciology* 49(164), 91–101.

Schulson EM (1990) The brittle compressive fracture of ice. *Acta Metallurgica Et Materialia* 38(10), 1963–1976. doi: 10.1016/0956-7151(90)90308-4

Schulson EM (2001) Brittle failure of ice. *Engineering Fracture Mechanics* 68(17), 1839–1887. doi: 10.1016/S0013-7944(01)00037-6

Schweizer J (1998) Laboratory experiments on shear failure of snow. *Annals of Glaciology* 26, 97–102. doi: 10.3189/1998AoG26-1-97-102

Shoop S, Kestler K and Haehnel R (2006) Finite element modeling of tires on snow. *Tire Science and Technology* 34(1), 2–37. doi: 10.2346/1.2169827

Smilauer V (2010) Cohesive particle model using the discrete element method on the Yade platform. PhD Thesis, Université de Grenoble; Czech Technical University in Prague.

Szabo D and Schneebeli M (2007) Subsecond sintering of ice. *Applied Physics Letters* 90(15), 151916. doi: 10.1063/1.2721391

Takei I and Maeno N (2004) Mechanical vibration responses of snow samples near the melting temperature. *Annals of Glaciology* 38, 130–134. doi: 10.3189/172756404781815022

Trottet B and 6 others (2022) Transition from sub-Rayleigh anticrack to supershear crack propagation in snow avalanches. *Nature Physics* 18(9), 1094–1098. doi: 10.1038/s41567-022-01662-4

Veylon G, Nicot F, Zhu H and Darve F (2018) Microstructure incidence on the bifurcation domain topology in granular materials. *Journal of Engineering Mechanics* 144(6), 04018031. doi: 10.1061/(ASCE)EM.1943-7889.0001446

Wautier A and 6 others (2021) Multiscale modelling of granular materials in boundary value problems accounting for mesoscale mechanisms. *Computers and Geotechnics* 134, 104143. doi: 10.1016/j.compgeo.2021. 104143

Wautier A, Bonelli S and Nicot F (2017) Scale separation between grain detachment and grain transport in granular media subjected to an internal flow. *Granular Matter* 19(2). doi: 10.1007/s10035-017-0706-9

Wautier A, Geindreau C and Flin F (2015) Linking snow microstructure to its macroscopic elastic stiffness tensor: A numerical homogenization method and its application to 3-D images from X-ray tomography. *Geophysical Research Letters* 42(19), 8031–8041. doi: 10.1002/2015GL065227

Willibald C, Scheuber S, Löwe H, Dual J and Schneebeli M (2019) Ice spheres as model snow: Tumbling, sintering, and mechanical tests. *Frontiers in Earth Science* 7, 229. doi: 10.3389/feart.2019.00229

Xiong H and 7 others (2021) A novel multi-scale large deformation approach for modelling of granular collapse. *Acta Geotechnica* 16(8), 2371–2388. doi: 10.1007/s11440-020-01113-5

Xiong H, Nicot F and Yin ZY (2017) A three-dimensional micromechanically based model. *International Journal for Numerical and Analytical Methods in Geomechanics* 41(17), 1669–1686. doi: 10.1002/nag.2692

Xiong H, Yin Z-Y and Nicot F (2019) A multiscale work-analysis approach for geotechnical structures. *International Journal for Numerical and Analytical Methods in Geomechanics* 43(6), 1230–1250. doi: 10.1002/nag. 2893

Zhu H, Nicot F and Darve F (2016a) Meso-structure organization in two-dimensional granular materials along biaxial loading path. *International Journal of Solids and Structures* 96. doi: 10.1016/j.ijsolstr.2016. 06.025

Zhu H, Veylon G, Nicot F and Darve F (2016) On the mechanics of meso-scale structures in two-dimensional granular materials. *European Journal of Environmental and Civil Engineering* 21(7–8), 912–935. doi: 10.1080/ 19648189.2016.1229229

Zhu H-X, Nicot F and Darve F (2016b) Meso-structure evolution in a 2D granular material during biaxial loading. *Granular Matter* 18, 1–12.

# Appendices

## A Expression of the tangential force in a bonded contact between ice grain

The total shear strain $\epsilon_{b,s}^{tot}$ in the bond can be split into an elastic part $\epsilon_{b,s}^e$ and a viscoplastic part $\epsilon_{b,s}^{vp}$:

$$ \epsilon_{b,s}^{tot} = \epsilon_{b,s}^{vp} + \epsilon_{b,s}^{e} \quad (39) $$

The viscoplastic deformation accounts for the creep behaviour of ice. Then, we assume that the superposition principle is applicable in Glen's law, as it has been done in Kabore and others (2021). This assumption enables us to express the normal and shear viscous strain independently. The normal creep rate is given as a function of the normal stress with the standard expression of Eqn (19) (Barnes and others, 1971). Then, for the viscous shear strain $\dot{\epsilon}_{b,s}^{vp}$, we can write:

$$ \dot{\epsilon}_{b,n}^{vp} = A(T_{\text{emp}}) \sigma_{b,s}^{m} \quad (40) $$

where $A(T_{\text{emp}})$ and $m$ are the same constant as in Eqn (19) and $\sigma_{b,s}$ is the shear stress in the bond.

In addition, the elastic shear deformation is linearly related to the shear stress

$$ \sigma_{b,s} = G \epsilon_{b,s}^{e} \quad (41) $$

with G the shear modulus of ice.

This gives the total viscous strain as a function of the shear stress:

$$ \epsilon_{b,s}^{tot} = \frac{1}{G} \sigma_{b,s} + A(T_{\text{emp}}) \sigma_{b,s}^{m} \quad (42) $$

Going back to the beam representation of the ice bond in the H-cell, we can express the shear strain in the bond as:

$$ \epsilon_{b,s}^{tot} = \frac{(\alpha - \alpha_0) L_b}{L_{b,0}} \quad (43) $$

with $(L_b, \alpha)$ and $(L_{b,0}, \alpha_0)$ being the bond length and the H-cell opening angle, respectively, at the current time and at the beginning of the loading (Kabore and others, 2021) (Peters and others, 2021).

And the shear stress can be related to the tangential contact force $T^b$ as follows:

$$ \sigma_{b,s} = \frac{T^b}{S_b} = \frac{T^b}{\pi r_b^2} \quad (44) $$

where $r_b$ is the cylindrical bond radius at time $t$.

Eventually, due to ice incompressibility in the bonds, the volume of an ice bond is supposed to be constant, so we get:

$$ \pi r_b^2 L_b = \pi r_{b,0}^2 L_{b,0} \quad (45) $$

By combining Eqns (39)–(45), the contact law relating the incremental normal contact force $\delta N^b$ and the incremental intergranular distance $\delta d$ reads:

$$ \delta T^b = - \frac{T \delta d}{L_b} + \pi r_b^2 G \frac{L_b \delta \alpha}{L_b (t_0)} - \pi r_b^2 G A(T_{\text{emp}}) \sigma_{b,t}^m \delta t \quad (46) $$

## B Contact stiffness after bond failure

After bond failure, we consider that the contact force is controlled by the overlapping of the two grains at contact. As the ice grains remained deformable, the contact stiffness is strongly dependent on the contact area between the two grains. The normal contact force $F_{n,i}$ is linked to the normal stiffness $k_n$ by the relationship:

$$ F_{n,i} = k_n (2r_g - d_i) \quad (47) $$

with $r_g$ being the grain radius and $d_i$ the intergranular distance.
The contact force can also be written as a function of the normal contact stress $\sigma_{n,i}$:

$$F_{n,i} = \sigma_{n,i} * (\pi r_c^2) \quad (48)$$

with $r_c = \sqrt{r_g^2 - d_i^2}/4$ being the radius of the contact between the grain. The normal stress can be expressed in terms of normal strain $\epsilon_n$:

$$\sigma_{n,i} = E_{\text{ice}} \epsilon_n = E_{\text{ice}} \frac{2r_g - d_i}{2r_g} \quad (49)$$

So, we finally get the following expression for the normal stiffness:

$$k_n = \frac{E_{\text{ice}}}{2r_g} * \pi r_c^2 = \frac{E_{\text{ice}} \pi}{2} \left( r_g - \frac{d_i^2}{4r_g} \right) \quad (50)$$

## C Creation of new contacts in the 3-D H-model H-cell

In the general case, the coordination number in a H-cell of the 3-D H-model is $Z = 2.4$. However, depending on the loading path, the geometry of the H-cell can be deformed to a point where additional contacts are created. We can distinguish two cases: the new contact is between extremal grains 1 and 4 ($Z = 2.6$) or between the two hexagons of the H-cell (contact between grains 2 and 7, $Z = 4$).

### C-1 Creation of contact between grains 1 and 4

The contact created between grains 1 and 4 generates an additional normal contact force $N_5$ to be added in the system and that contributes to the stresses in the H-cell. The value of the contact force depends on the intergranular distance $d_5$ between grains 1 and 4, as defined in Fig. 17 and reads:

$$N_5 = - \max(k_n(2r_g - d_5); 0) \quad (51)$$

with $k_n$ depending on the contact law used in the model. Adding this force does not impact the equilibrium on grain 2 and consequently does not change the way to obtain the variations of the parameters of the H-cell geometry. Only the expression of the local axial stress must be modified as follows:

$$\tilde{\sigma}_{11} = \sigma_{11}^{\text{other contacts}} + \frac{N_5 d_5}{V_{\text{meso}}} \quad (52)$$

with $\tilde{\sigma}_{11}^{\text{other contacts}}$ the contribution of the other contact forces on the first component of the local stress tensor.

At some point, the opening angle may reach 90°, which means that its value will not be able to increase anymore. Additional external forces must be added on grain 2 to prevent any further variation of the opening angle (see Fig. 18). In this case, the evolution of the geometry of the H-cell is directly given by $\delta L_1$, $\delta L_2$ and $\delta L_3$:

$$\left\{ \begin{array}{l} \delta d_1 = \delta L_2/2 \\ \delta d_2 = \delta d_5 = \delta L_1 \\ \delta \alpha = 0 \end{array} \right. . \quad (53)$$

<figure><img src="image_19.png" /><figcaption>**Figure 17.** Description of the creation of contact between grains 1 and 4, with $\alpha < 90°$: (a) scheme of one hexagon of the H-cell; (b) force equilibrium on grain 1.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>**Figure 18.** Description of the creation of contact between grains 1 and 4, with $\alpha = 90°$: (a) scheme of one hexagon of the H-cell; (b) force equilibrium on grain 2.</figcaption></figure>

### C-2 Creation of contact between the two hexagons of the same H-cell

When the lateral strains become high enough, new contacts between the two hexagons of the 3-D H-cell can appear. Considering the possibility that $d_2(t) \neq d_4(t)$, the direction of the contact can be characterized by two Eulerian angles $\phi_3$ and $\theta_3$ as defined in Fig. 19. $\phi_3$ is the angle between the x-axis and the branch vector $\bar{d}_6$ between the centres of grains 2 and 7 while $\theta_3$ is the angle between the z-axis and the orthogonal projection of the branch vector $\bar{d}_6$ on the (Oyz) plane. These angles can be written:

$$\left\{ \begin{array}{l} \phi_3 = \arctan \frac{\sqrt{d_2^2 + d_3^2}}{d_3 \sqrt{d_2^2 - d_3^2}} \\ \theta_3 = \arctan \frac{d_2 \sqrt{d_2^2 - d_3^2}}{z_7} \end{array} \right. \quad (54)$$

with $y_2$ the y-coordinate of grain 2 and $z_7$ the z-coordinate of grain 7.

This new contact induces the addition of two new contact forces $N_6$ and $T_6$ between grains 2 and 7 (see Fig. 20). The value of these forces depends on the length of the branch vector $\bar{d}_6$ given by the following expression:

$$\bar{d}_6 = \begin{pmatrix} d_4/2 - d_2/2 \\ d_1 \sin \alpha_1 \\ d_3 \sin \alpha_2 \end{pmatrix}. \quad (55)$$

From this it follows that:

$$d_6 = \left( \left( d_4/2 - d_2/2 \right)^2 + (d_1 \sin \alpha_1)^2 + (d_3 \sin \alpha_2)^2 \right)^{1/2} . \quad (56)$$

The contact force $\bar{N}_6$ and the vector $\bar{d}_6$ being collinear, we can write:

$$\overline{N_6(t)} = k_n(2r_g - d_6(t)) \frac{\overline{d_6(t)}}{d_6(t)} = \frac{k_n(2r_g - d_6(t))}{d_6(t)} \begin{pmatrix} d_4/2 - d_2/2 \\ d_1 \sin \alpha_1 \\ d_3 \sin \alpha_2 \end{pmatrix} \quad (57)$$

with $k_n$ depending on the contact law used in the model. In the case of a frictional law, the magnitude of the tangential contact force $\bar{T}_6$ reads:

$$|T_6(t)| = \min |T_6(t - \delta t) + k_n d_6 \delta \alpha_3|, \quad |N_6(t) \tan \phi_g| \quad (58)$$

where $\delta \alpha_3(t)$ is the angle between $\overline{d_6}(t - \delta t)$ and $\overline{d_6}(t)$, defined by:

$$\delta \alpha_3 = \arccos \frac{\langle \overline{d_6}(t) | \overline{d_6}(t - \delta t) \rangle}{\| \overline{d_6}(t) \| \cdot \| \overline{d_6}(t - \delta t) \|} \quad (59)$$

The tangential contact force $\bar{T}_6(t)$ is orthogonal to the normal contact force $\overline{N_6}(t)$ and coplanar with $\overline{N_6}(t)$ and $\overline{N_6}(t - \delta t)$, which means that it verifies:

$$\left\{ \begin{array}{l} \langle \overline{T}_6(t) | \overline{N}_6(t) \rangle = 0 \\ \det(\overline{T}_6, \overline{N}_6(t), \overline{d}_6(t - \delta t)) = 0 \end{array} \right. . \quad (60)$$

So $\overline{T}_6$ can be written:

$$\overline{T}_6 = \begin{pmatrix} T_{6x} \\ T_{6y} \\ T_{6z} \end{pmatrix} \quad (61)$$
<figure><img src="image_21.png" /></figure>

<figure><img src="image_22.png" /><figcaption>**Figure 19.** Description of the creation of contact between grains 2 and 7, with $d_2 \neq d_4$: (a) scheme of the H-cell and definition of the $\phi_3$ angle; (b) projection of the new contact in the plan (Oyz) and definition of the angle $\theta_3$.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>**Figure 20.** Force equilibrium on grain 2: (a) in the plan (Oxy); (b) in the plan (Ozy).</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>**Figure 21.** Description of the loss of contacts: (a) loss of the inclined contact between grains 1 and 2; (b) and (c) loss of the axial contact between grains 2 and 3 with (b) the description of the equilibrium on the two half cells and (c) the description of equilibrium on grains 1 and 2.</figcaption></figure>

with

$$ \left\{ \begin{aligned} T_{6x} &= a_{tx} T_{6y} \\ T_{6y} &= \frac{|\overline{T_6(t)}|}{\sqrt{1+a_{tx}^2+a_{ty}^2/a_{tz}^2}} \\ T_{6z} &= \frac{a_{fy}}{a_{tz}} T_{6y} \end{aligned} \right. \qquad (62) $$

and

$$ \left\{ \begin{aligned} a_{tx} &= \frac{-d_{6y}(t)}{d_{6x}(t)} - \frac{d_{6z}(t)}{d_{6x}(t)} \frac{a_{ty}}{a_{tz}} T_{6y}(t) \\ a_{fy} &= \frac{(d_{6y}(t)d_{6z}(t-\delta t)-d_{6z}(t)d_{6y}(t-\delta t))}{d_{6x}(t)} + (d_{6x}(t)d_{6z}(t-\delta t)-d_{6z}(t)) \\ a_{tz} &= (d_{6x}(t)d_{6y}(t-\delta t)-d_{6y}(t)d_{6x}(t-\delta t)) - \frac{d_{6z}(t)}{d_{6x}(t)} \\ & \times (d_{6y}(t)d_{6z}(t-\delta t)-d_{6z}(t)d_{6y}(t-\delta t)) \end{aligned} \right. \qquad (63) $$

For the sake of simplicity and to keep the scheme of the H-model explicit, we do not account for the influence of $\overline{N_6}$ and $\overline{T_6}$ on the force equilibrium on grain 2, which means that the new forces only contribute to the expression of the stresses in the H-cell as follows:

$$ \left\{ \begin{aligned} \tilde{\sigma}_{11} &= \sigma_{11}^{\text{other contacts}} + \frac{8N_{6x}d_{6x}-8T_{6x}d_{6x}}{V_{\text{meso}}} \\ \tilde{\sigma}_{22} &= \sigma_{22}^{\text{other contacts}} + \frac{8N_{6y}d_{6y}-8T_{6y}d_{6y}}{V_{\text{meso}}} \\ \tilde{\sigma}_{33} &= \sigma_{33}^{\text{other contacts}} + \frac{8N_{6z}d_{6z}-8T_{6z}d_{6z}}{V_{\text{meso}}} \end{aligned} \right. \qquad (64) $$

## D Contact loss in the 3-D H-cell

Depending on the load path of a mesostructure, some contacts between grains in a H-cell can be lost. In the case of an unbonded contact, it means that

$d_1 > 2r_g$. In the case of a bonded contact, it occurs when the bond fails. We can distinguish two situations, depending on which contact is lost (Fig. 21). For the sake of simplicity, we consider in the following only one hexagon independently of the second one.

If the inclined contact (namely the one between grains 1 and 2) is lost (see Fig. 21a), the contact forces $N_1$ and $T_1$ are null, so the force equilibrium of grains 1 and 2 imposes that the external forces $F_1$ and $F_2$ are null too. The momentum equilibrium on grain 2 also imposes $T_1 = G_2 = 0$, which means that no force can equilibrate the normal contact force $N_2$ between grains 2 and 3. In that case, no force transmission is possible in the hexagon, so all the stresses are null. In this work, we have chosen to let the deformation of the cell continue until the contact is restored. We fix $d_2 = 2r_g$, and the geometrical compatibility equations (4) and (5) enables to determine the variation of $d_1$ and $\alpha$.

If the axial contact between grains 2 and 3 is lost, the normal contact force $N_2$ is null. In this case, one hexagon is separated into two triangles (see Fig. 21b, c). The equilibrium equation of forces and momentum on grain 2 leads to:

$$ \begin{gathered} G_2 = T_1 \\ N_1 \cos \alpha - T_1 (1 + \sin \alpha) = 0 \end{gathered} \qquad (65) $$

where $N_1$ and $T_1$ are of opposite signs. This means that the system to solve in this case becomes:
<figure><img src="image_25.png" /><figcaption>**Figure 22.** Description of the connectivity the grain in a cell with a loss of contact between grains 2 and 3 (grey cell), throughout the neighbouring cells.</figcaption></figure>

$$ \begin{pmatrix} 2 \cos \alpha_1 & 1 & -2d_1 \sin \alpha_1 \\ 2 \sin \alpha_1 & 0 & 2d_1 \cos \alpha_1 \\ A_1 & 0 & C_1 \end{pmatrix} \begin{pmatrix} \delta d_1 \\ \delta d_2 \\ \delta \alpha_1 \end{pmatrix} = \begin{pmatrix} \delta L_1 \\ \delta L_2 \\ \lambda_1 \end{pmatrix} \quad (66) $$

with $A_1$, $C_1$ and $\lambda_1$ as defined in Eqn (12), depending on the contact law.

We can calculate $F_1^a$ from the equilibrium on the forces on grain 1 projected on the axis of the cell:

$$ F_1^a = 2N_1 \cos \alpha - 2T_1 \sin \alpha \quad (67) $$

From Eqn (65), we deduce that:

$$ F_1^a = 2G_2 \quad (68) $$

which proves that the two half-cells are in equilibrium.

This external tangential force $G_2$ represents the action of the grains outside the cell on the grain within the cell, indicating that connectivity is established through contact with neighbouring cells. As illustrated in Fig. 22, the connectivity of grains 1 and 4 in the grey cell is ensured by their contact with the blue and orange cells.