# Revisiting snow settlement with microstructural knowledge

Louis Védrine, Pascal Hagenmuller

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

## To cite this version:

Louis Védrine, Pascal Hagenmuller. Revisiting snow settlement with microstructural knowledge. The Cryosphere, 2026, 20 (2), pp.1257-1277. 10.5194/tc-20-1257-2026. hal-05526610

HAL Id: hal-05526610

https://hal.science/hal-05526610v1

Submitted on 25 Feb 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished 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
<figure><img src="image_3.png" /></figure>

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

# Revisiting snow settlement with microstructural knowledge

Louis Védrine and Pascal Hagenmuller

Météo-France, CNRS, Univ. Grenoble Alpes, Univ. Toulouse, CNRM, Centre d'Études de la Neige, 38000 Grenoble, France

Correspondence: Louis Védrine (louis.vedrine@meteo.fr)

Received: 27 August 2025 – Discussion started: 17 September 2025

Revised: 12 December 2025 – Accepted: 28 January 2026 – Published: 18 February 2026

**Abstract.** Snow settlement under gravity is primarily driven by the slow creep of its ice matrix, which exhibits a viscoplastic behaviour. Knowledge of the viscoplastic properties of snow is thus crucial for understanding and predicting the snowpack seasonal and perennial evolution. However, different approaches have yielded disparate constitutive viscoplastic laws. Field experiments and snowpack models typically describe snow settlement with a linear model and an apparent compaction viscosity. Dedicated laboratory experiments exhibited non-linear relationships between stress and strain rate, with a stress exponent ranging from 1.8 to 4. Microstructure-based simulations showed that the viscoplastic behaviour likely results from the interaction of glides on various intra-crystalline slip systems within ice crystals, yielding an exponent between 2 and 3. The paper aims to reconcile these approaches. To do so, we conducted microstructure-based simulations on 37 three-dimensional snow images and established that the viscoplastic behaviour follows a power-law relation, with a stress exponent almost constant around 2.15, and a reference stress that depends mostly on the solid fraction. Analysing a dataset from previous viscoplastic tests (178 points) revealed that applying the stress exponent from the simulations significantly reduces variability in the reference stress between independent studies and led to a simplified constitutive relation. Lastly, we showed that the linear settlement laws of snowpack models, such as Crocus, align with the proposed constitutive relation under natural loading conditions typically encountered on alpine sites, due to correlations between stress and density. However, considerable differences emerge under “non-standard” scenarios, such as elevated loads on light snow or reduced loads on dense snow, where our model demonstrates superior qualitative performance.

## 1 Introduction

Under its weight, snow naturally settles from an initial density of approximately 100 kg m⁻³ just after deposition on the ground to about 500 kg m⁻³ at the end of the winter season, and up to the ice density for perennial snow on glaciers and ice sheets. Snow settlement is a key process of the snow-pack evolution (Brun et al., 1992; Lehning et al., 2002; Simpson et al., 2021). Predicting snow density is thus essential for many applications, such as avalanche forecasting (Morin et al., 2020), hydrology (Barnett et al., 2008), surface mass balance in polar regions (Morris and Wingham, 2014), or palaeoclimatology (Barnola et al., 1987).

Dry snow is a two-phase material made of sintered ice crystals and air (Roscoat et al., 2011). Its settlement is primarily driven by the slow viscoplastic deformation of its ice skeleton, leading to a progressive reduction in pore volume (e.g., Yosida, 1955). Indeed, snow is an open-pore material and the viscous flow of air during deformation can be neglected at low strain rates (e.g., Bergen, 1968; Colbeck, 1989). Snow settlement thus depends solely on the spatial arrangement of the ice crystals and the deformation mechanisms they undergo. Two main mechanisms are generally considered for ice in snow: intra-crystalline deformation, also known as dislocation glide (Duval et al., 1983; Louchet, 2004; Weikusat et al., 2017), and inter-crystalline deformation, also known as grain boundary sliding (Alley, 1987; Lundin et al., 2017). At macroscale, in uniaxial compression, the viscoplastic behaviour of snow, or more widely of porous polycrystals, is generally modelled by a power-law relating stress σ to strain rate ˙ε (Glen, 1955):

$$ \dot{\epsilon} = \dot{\epsilon}_0 \left( \frac{\sigma}{\sigma_0} \right)^n \qquad (1) $$
with $\dot{\epsilon}_0 = 1 \text{ s}^{-1}$, n the stress exponent and $\sigma_0$ the reference stress. This model is extremely sensitive to n, which quantifies the dependence of the viscous deformation on the applied stress.

To interpret experimental compression tests on snow, the concept of compaction viscosity has often been employed (Bader, 1954, 1962; Mellor, 1974; Navarre, 1975; Gray and Morland, 1995; Kominami et al., 1998; Sturm and Holmgren, 1998; Scapozza and Bartelt, 2003a). The compaction viscosity $\eta$ linearly relates the two directly-measurable macroscopic quantities, i.e., stress $\sigma$ and strain rate $\dot{\epsilon}$:

$$ \dot{\epsilon} = \frac{\sigma}{\eta}. \qquad (2) $$

This model implicitly assumes that $n = 1$ in Eq. (1) and corresponds to a Newtonian compressible viscous fluid. Several authors have encapsulated into $\eta$ various processes presumably contributing to snow creep, such as gravitational settling, grain rearrangement, water infiltration, and refreezing (Shinojima, 1967; Kojima, 1967; Bergen, 1968; Keeler, 1969). Rheological models, for instance, the Maxwell, Kelvin, or Burger models (a combination of linear springs and dampers in parallel or series) are also used to reproduce the permanent viscous flow at low strain rates (Bader, 1962; Mellor, 1974; McClung and Larsen, 1989; Bartelt and Christen, 1999; Chandel et al., 2007; Gorynina et al., 2024; Bahaloo et al., 2024; Huo et al., 2024). Nevertheless, these models apply to a wide range of conditions only if $\eta$ does not depend on $\sigma$ and $\dot{\epsilon}$, which remains questionable. Consequently, non-linear extensions of these viscoelastic models have been proposed (Desrues et al., 1980; Stoffel and Bartelt, 2003; O'Connor and Haehnel, 2020). Besides, in line with this historical data and models, the settlement laws currently implemented in snowpack models, such as Crocus (Brun et al., 1992; Vionnet et al., 2012) and Snowpack (Lehning et al., 2002; Simson et al., 2021), are also linear. In these models, the compaction viscosity is parametrized as a function of temperature and snow microstructure (primarily its density) and has been calibrated and evaluated under specific conditions (Ganju et al., 1999). Again, it is questionable whether the implemented models apply to conditions that deviate from the calibration conditions, typically a natural snowpack evolving under alpine conditions. For instance, Royer et al. (2021), Woolley et al. (2024) have tuned the standard settlement laws to represent the evolution of snowpack density under arctic conditions, but they lack a physical basis and a wealth of calibration data.

Different laboratory experiments have focused on characterizing n on snow. Narita (1980, 1984) reported a stress exponent of $n = 4 \pm 1$ which appeared independent of temperature. Scapozza and Bartelt (2003a, b) conducted strain rate-controlled triaxial experiments on alpine snow and identified two distinct regimes: regime I with an exponent of approximately $n \approx 1.8$ for low-density snow, and regime II with $n \approx 3.5$ for higher density snow, with a transition

density depending on temperature. In line with this finding, Schleef et al. (2014a) reported $n \approx 2$ for compression experiments on low-density samples. Delmas (2013) observed a sensitivity of n to snow types: $n = 4.1$ for faceted grains compared to $n = 3.2$ for rounded grains, at $-15^{\circ}$. Recently, Sundu et al. (2024b) performed cyclic loadings on snow samples and derived n from the loading and relaxation phases. They found a transition of n as a function of the geometrical grain size: from $n \approx 1.9$ for small grains to $n \approx 4.4$ for large grains. They attributed this transition to a transition from grain boundary sliding to dislocation creep at a critical grain size. In general, the reported scatter of n between 1.8 and 4 was often attributed to the driving mechanisms at the micro-scale, whose activation may depend on density, optical diameter, or snow type. However, experimental evidence for the reasons behind this sensitivity is still lacking (Ignat and Frost, 1987; Meyssonnier et al., 2009). Moreover, repeated experiments under perfectly controlled conditions (e.g., changing the applied stress without altering any parameter that affects the reference stress $\sigma_0$) are seldom conducted and are challenging on snow.

An approach to understanding the mechanisms of snow viscoplasticity and predicting the overall behaviour relies on numerical experiments. The key idea is to model snow based on its porous and crystalline microstructure and to incorporate the local deformation mechanisms. The numerical solver then simulates the macroscopic behaviour of the sample. Thanks to higher computing power and easier access to microstructural data, various studies have followed this approach in recent years for snow viscoplasticity, with different assumptions regarding the geometry of the ice skeleton and the driving mechanisms (Johnson and Hopkins, 2005; Theile et al., 2011; Wautier et al., 2017; Védrine et al., 2024, 2025). Johnson and Hopkins (2005) modelled snow creep with a snow microstructure described with a set of cylinders with hemispherical ends and inter-granular viscoplastic contact law. Theile et al. (2011) and Wautier et al. (2017) employed a more accurate representation of the microstructure based on tomographic data, and considered only intracrys talline deformation. These studies simulated an exponent $n = 3$ and achieved reasonable agreement with experimental data. However, the validation data were either too far from the numerical conditions to perform a 1:1 evaluation (Wautier et al., 2017) or too scarce to validate the model assumptions definitely (Theile et al., 2011). In addition, a number of their modelling assumptions are now questionable a priori. Theile et al. (2011) modelled the snow microstructure as a collection of interconnected finite element beams, which inadequately captured the deformation obstacles present at the crystal boundaries. Conversely, Wautier et al. (2017) treated the ice in snow as a homogeneous isotropic polycrystalline material, neglecting the fact that ice in snow is a collection of ice crystals surrounded by air. Indeed, Védrine et al. (2024) recently demonstrated that ice in snow cannot be modelled as an isotropic material, and that the crystalline nature of
ice must be explicitly considered. In general, Védrine et al. (2025) showed that, for porous polycrystals like snow, the stress exponent is sensitive to density and intracrystalline surface area. They showed that the homogenized stress exponent varies between the exponent of the weakest (basal) slip system ($n = 2$, Duval et al., 1983; Chevy, 2008) and the stress exponent of dense polycrystalline ice ($n \sim 3$), depending on the frustration of the crystals. Here, “frustration” refers to the mechanical incompatibility between neighbouring crystals that constrains soft deformation mechanisms.

In summary, three different approaches have been used to understand and model the viscoplastic behaviour of snow, which is essential to predict its settlement:

- Field and laboratory experiments focusing on conditions typically encountered by natural snowpacks in alpine regions have described snow settlement by a linear model ($n = 1$) and an apparent compaction viscosity $\eta$ (Eq. 2). This law is implemented in detailed snowpack models.

- Laboratory experiments focusing on the generic material properties of snow reported a non-linear relationship between stress and strain rate (Eq. 1) with a stress exponent $n$ between 1.8 and 4.

- Microstructure-based simulations suggested that the macroscopic viscoplastic behaviour results from the complex interplay between glides on different ice intracrystalline slip systems, which yields an exponent $n$ between 2 and 3.

Our paper aims to show that, despite these apparently contradictory results, these three approaches can be reconciled. Our starting point is the simulations of the effective viscoplastic behaviour of numerous snow microstructures with the state-of-the-art 3D model by Védrine et al. (2025). Then, we re-interpret previously published experimental data by searching for the model with the minimal scatter. We show that the model ends up very close to the one derived from microstructure-based simulations. Eventually, we analyse the difference between the latter model and the settlement law implemented in Crocus and show that they coincide under the conditions typically encountered by natural snowpacks in the Alps.

# 2 Material and Methods

## 2.1 Microstructure-based simulations

We conducted numerical full-field simulations on three-dimensional porous microstructures, using a crystal plasticity model to describe the local deformation mechanisms of ice. These simulations are then homogenized to derive the stress exponent $n$ and reference stress $\sigma_0$ (Eq. 1) of a given snow microstructure.

### 2.1.1 Input snow microstructures

The input of our microstructure-based viscoplastic model is a 3D image of the snow microstructure partitioned into individual crystals. We used 37 different X-ray tomographic images acquired from previously published studies (Hagenmuller et al., 2016; Peinke et al., 2020; Hagenmuller and Car magnola, 2020; Fourteau et al., 2022; Bernard et al., 2023) (Table 1). These binary ice-air images encompass the majority of snow types encountered in seasonal snowpacks and span a wide range of density and specific surface area (SSA), i.e., the ice-air surface area per unit snow mass (Fig. 1). X-ray attenuation images do not distinguish between ice crystals. We thus numerically partitioned the ice skeleton into individual grains based on geometrical criteria, following the algorithm of Hagenmuller et al. (2014b) and Peinke et al. (2020). We assumed that the detected grains correspond to the ice crystals, which is a reasonable assumption for most types of snow according to Arnaud et al. (1998). Then, each crystal was associated with a full crystallographic orientation (c-axis and a-axes) randomly sampled in an isotropic distribution. Note that Védrine et al. (2024) showed that the exact random draw does not significantly affect the overall mechanical behaviour. To reduce computation time, simulations were not performed on the entire scanned volume at the nominal image resolution. The resolution was downscaled while ensuring that the equivalent optical radius, $3/(SSA\rho_{ice})$, remains between 6 and 12 voxels. A cubic sub-volume of $250 \times 250 \times 250$ voxels was extracted from the image. This volume contains at least 20 crystals along each edge, which ensures the sample is representative (Védrine et al., 2024). Some input 3D images are shown as 2D slices in Fig. 2. All details on these images are given in the supplementary material and associated publications.

### 2.1.2 Elasto-viscoplastic simulations

We performed elasto-viscoplastic simulations using the model developed by Védrine et al. (2024, 2025). We recall here its main assumptions. Details can be found in the Appendix A and associated references.

Air was modelled as an infinitely soft elastic medium. An elasto-viscoplastic law modelled the mechanical behaviour of each ice crystal. The elasticity tensor was assumed isotropic transverse in the crystal frame and defined according to Gammon et al. (1983). The viscoplastic part follows the crystal plasticity law proposed by Castelnau et al. (2008), Lebensohn et al. (2009), and Suquet et al. (2012). Ice crystals deform through slip on three soft basal systems, three hard prismatic systems, and six hard pyramidal systems (Montagnat et al., 2014; Weikusat et al., 2017). Each slip system ($k$) is governed by a Norton-type flow law, characterised by a prefactor, the reference shear rate $\dot{\gamma}_0^{(k)}$, a stress exponent $n^{(k)}$, and a critical shear stress $\tau_0^{(k)}$. Different parametrizations of the slip systems have been proposed in the literature for dif-
<table><caption>**Table 1.** Sample origin and reference papers, type and number of samples and voxels size (range) of images used in this study.</caption><thead><tr><th>Origin and reference</th><th>Type</th><th>Numbers of scans</th><th>Resolution (μm)</th></tr></thead><tbody><tr><td>Hagenmuller et al. (2016)</td><td>RG, DH</td><td>8</td><td>[10,18]</td></tr><tr><td>Peinke et al. (2020)</td><td>PP, RG, RGlr, DH</td><td>9</td><td>15</td></tr><tr><td>Hagenmuller and Carmagnola (2020)</td><td>MM</td><td>2</td><td>15</td></tr><tr><td>Fourteau et al. (2022)</td><td>DH, RG, MF, FC</td><td>8</td><td>15</td></tr><tr><td>Bernard et al. (2023)</td><td>DF</td><td>10</td><td>8</td></tr></tbody></table>

The type of snow is defined according to Fierz et al. (2009). The numbers shown in brackets, e.g., [a, b], mean that the parameter ranges between, e.g., a and b. Details can be found in the associated references.

<figure><img src="image_5.png" /><figcaption>**Figure 1.** Characteristics of the 3D microstructures used in this study. The snow types are described and coloured according to the International Classification for Seasonal Snow on the Ground (Fierz et al., 2009): Precipitation Particles (PP), Decomposing and Fragmented Precipitation Particles (DF), Rounded Grains (RG), Large Rounded Grains (RGr), Faceted Crystals (FC), Depth Hoar (DH), Machine Made (MM). The microstructures enclosed by a black circle correspond to those shown in Fig. 2.</figcaption></figure>

ferent purposes: to investigate the viscoplastic development of texture (Castelnau et al., 1996), to model the viscoplastic behaviour and heterogeneous intracrystalline deformation of columnar ice polycrystals (Lebensohn et al., 2009), or to simulate the transient regime (Suquet et al., 2012; Castelnau et al., 2008). However, as shown by Védrine et al. (2025), a key point for modelling snow behaviour is to set the exponent for the basal system equal to 2, as extensively experimentally observed (Duval et al., 1983; Schulson and Duval, 2009). We therefore used the stationary parametrization of Suquet et al. (2012), calibrated to reproduce the viscoplastic behaviour of dense polycrystalline ice ($\sigma_0(\Phi = 1) = 272$ MPa and $n_{suquet} = 2.94$) as proposed by Budd and Jacka (1989) and Castelnau et al. (1996) at $T_0 = 263.15$ K (for more details see Appendix A). To account for temperature effects, we used the Arrhenius relation proposed by Mellor and Testa

(1969a), i.e., $\dot{\gamma}_0^{(k)}(T) = \dot{\gamma}_0^{(k)}(T_0)A(T)$ with

$$A(T) = \exp \left[ - \frac{Q}{R} \left( \frac{1}{T} - \frac{1}{T_0} \right) \right], \quad (3)$$

where $Q = 69.1$ kJmol$^{-1}$ is the activation energy of polycrystalline ice (Mellor and Testa, 1969b), $R = 8.3$ Jmol$^{-1}$K$^{-1}$ is the universal gas constant, $T$ is the absolute temperature (in kelvin), and $T_0 = 263$ K denotes the reference temperature at which the parameters are defined.

For the boundary conditions, we prescribed vertical uniaxial compression with periodic boundary conditions (Fig. 3). We imposed a positive average vertical strain rate $\dot{\varepsilon}_{zz} > 0$ and zero average lateral stresses, $\sigma_{xx} = \sigma_{yy} = 0$, and zero average shear stresses, $\sigma_{xy} = \sigma_{yz} = \sigma_{xz} = 0$. Here, the macroscopic strain and stress tensors are defined as volume average on the sample: $\varepsilon = \frac{1}{|V|} \int_V \varepsilon(x)dV$ and $\sigma = \frac{1}{|V|} \int_V \sigma(x)dV$, respectively (Hill, 1963). These boundary conditions mimic the conditions of snowpack densification, where gravity is the primary driving force and the influence of confining pressure is negligible (Wautier et al., 2017). We measured the resulting average vertical stress $\sigma_{zz}$ from the simulations. By convention, as the problem is purely compressive, both the strain rate and stress are defined as positive.

The simulations were conducted using the Fast Fourier Transform (FFT)-based solver AMITEX_FFTP (Gelebart et al., 2020), using the small perturbation assumption. This solver employs a modified Green operator based on a hexahedral finite difference scheme (Willot, 2015), equivalent to using linear hexahedral finite elements with reduced integration (Schneider et al., 2017), and an Anderson convergence acceleration technique (Anderson, 1965), enabling the analysis of porosities (infinite contrast with the solid matrix) up to 90%. The crystal plasticity model was implemented using the MFront code generator (Helfer et al., 2015). The numerical integration of crystal plasticity is computationally expensive (approximately 12 h on the sample presented in Fig. 3 of size $250 \times 250 \times 250$ voxels, on a high-performance computer with 256 IntelXeon Gold 6130 cores). However, the solver benefits from the MPI implementation.
<figure><img src="image_6.png" /><figcaption>**Figure 2.** Example of the *c*-axis orientation of different samples: (a) PP, (b) RG, (c) DF, (d) RGlr, (e) DH, (f) MM. Each colour represents a specific *c*-axis orientation, as indicated by the accompanying colour wheel. The black background corresponds to the original air pores within the snow microstructure. Note that the simulation volumes are 3D. Here, only 2D slices are shown.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 3.** Numerical boundary conditions: uniaxial compression test with imposed strain rate.</figcaption></figure>

### 2.1.3 Homogenized model

The mechanical behaviour of heterogeneous materials, such as snow, can often be modelled by replacing the material with an equivalent homogeneous medium. This approximation,

also known as homogenization, requires that the characteristic length scale of the material microstructure is small relative to the macroscopic scale of interest (Dormieux and Bourgeois, 2002; Auriault et al., 2010). In the case of snow, this separation of scales is typically satisfied, which allows for the macroscopic mechanical behaviour to be deduced from the mesoscopic structure obtained through X-ray tomography. Previous studies have shown that snow samples on the order of a few millimeters can be considered representative elementary volumes (REVs) for modelling the mechanical behaviour of snow (Wautier et al., 2015; Srivastava et al., 2016; Védrine et al., 2024; Sundu et al., 2024a).

Figure 4 illustrates the uniaxial compression of the exemplary microstructure shown in Fig. 3, under strain rates ranging from $10^{-7}$ to $10^{-5}$ s$^{-1}$. We observe three different regimes of the stress-strain response (Fig. 4). First, the axial stress $\sigma_{zz}$ increases almost linearly with strain in an elastic regime. Indeed, at very low stress levels and short times, the viscoplastic strain is negligible. Then, $\sigma_{zz}$ deviates from the elasticity, and the viscoplasticity is progressively activated. After this transient regime, $\sigma_{zz}$ eventually reaches a constant value, the flow stress $\sigma_y$, where the microstructure flows perfectly (stationary creep). The flow stress depends on the microstructure and the applied strain rate. In our simulations, the permanent viscoplastic regime is reached after a typi-
<figure><img src="image_8.png" /><figcaption>Figure 4. Evolution of the vertical stress $\sigma_{zz}$ as the function of the strain $\epsilon_{zz}$ for different strain rates. Inset: Evolution of axial flow stress as a function of the imposed strain rate. The tested sample is composed of Decomposing and Fragmented snow particles (DF) with a solid fraction $\Phi = 0.48$, and is shown in Fig. 3.</figcaption></figure>

cal strain of 0.5 %. This is consistent with the steady (minimum) strain rate of polycrystalline ice, commonly measured around 1 % (Mellor and Cole, 1982; Jacka, 1984; Treverrow et al., 2012). Consequently, we conduct simulations up to 1 % strain to determine the flow stress $\sigma_Y$.

The evolution of flow stress with strain rate follows a power-law over a wide range of strain, encompassing typical strain rates for snow settlement (Fig. 4 inset). Consequently, we can model the constitutive behaviour using the power-law of Eq. (1) by fitting two parameters $n$ and $\sigma_0$. The sensitivity to temperature implemented on the microscopic slip systems, $A(T)$ (Eq. 3), directly scales out at the macroscopic scale. In other words, the microstructure-based simulations yield the following model:

$$ \dot{\epsilon} = \dot{\epsilon}_0(T_0)A(T)\left(\frac{\sigma}{\sigma_0}\right)^n \quad (4) $$

with $n$ and $\sigma_0$ two fitting parameters depending on the tested microstructure.

## 2.2 Revisiting previous viscoplastic tests

In this section, we describe how to compare previously published viscoplastic tests and gather them into a common framework.

To characterize the viscoplastic behaviour of snow, two types of mechanical tests are commonly employed: load-controlled tests (creep tests) (Mellor and Testa, 1969b; Desrues et al., 1980; Bartelt and Christen, 1999; Kirchner et al., 2001; Scapozza and Bartelt, 2003a, b; Moos et al., 2003) and strain rate controlled tests (Scapozza and Bartelt, 2003a; Bernard, 2023). We selected experimental data from such tests, where the initial density, temperature, stress, and strain rate were measured during the experiments. To focus exclusively on the viscoplastic behaviour, only tests with a strain rate lower than $4 \times 10^{-4} \text{ s}^{-1}$ were selected (Kirchner et al., 2001; Barraclough et al., 2017; Löwe et al., 2020; Bernard, 2023). Our dataset comprises 92 measurements obtained from strain rate-controlled tests and 86 measurements obtained from stress-controlled tests (Scapozza and Bartelt, 2003b; Theile et al., 2011; Schleef et al., 2014b; Bernard, 2023; Bernard et al., 2023). One measurement comprises the measured strain rate and stress $(\dot{\epsilon}, \sigma)$ and essential information such as initial ice volume fraction $\Phi$, temperature $T$ and loading conditions (load or strain-controlled). The data is summed up in Table 2 and fully available in Védrine and Hagenmuller (2025).

It is impossible to derive two parameters, namely the reference stress $\sigma_0$ and stress exponent $n$, from an individual measurement point $(\dot{\epsilon}, \sigma)$. Moreover, the data points cannot be directly compared since they were obtained at various temperatures in the range [-18, -2] °C. However, given a certain $n$ in the range [1, 4] and assuming that the Arrhenius temperature scaling applies (Eq. 3), we rescale the data at a constant temperature of -10 °C, and we can calculate $\sigma_0(n)$ as:

$$ \sigma_0(n) = \sigma \times \left( \frac{\dot{\epsilon}}{\dot{\epsilon}_0(T_0)A(T)} \right)^{-\frac{1}{n}} \quad (5) $$

If the power-law is well-formed, $\sigma_0$ and $n$ do not depend on the applied stress or strain rate but solely rely on the sample microstructure. Therefore, for a given microstructure, knowing the correct stress exponent $n$ is equivalent to obtaining the same reference stress regardless of the test conducted (loading level or test type). In other words, finding the appropriate constitutive law is equivalent to identifying the exponent $n$ yielding the smallest scatter (defined as the relative standard deviation) of $\sigma_0(n)$ for different experimental tests but similar microstructures.

In this study, we proposed a simplified but pragmatic application of this procedure by considering the solid fraction $\Phi$ as the sole characteristic of the microstructure. Indeed, the characterization of microstructures from data collected in the literature is primarily limited to density. Therefore, we assumed $\sigma_0 = \sigma_0(n, \Phi)$ and $n = n(\Phi)$. Consequently, the ex-
<table><caption>**Table 2.** Characteristic of the experimental data of snow creep used.</caption><thead><tr><th>References</th><th>Snow type</th><th>Temperature (°C)</th><th>Density (kg m<sup>-3</sup>)</th><th>Strain rate (s<sup>-1</sup>)</th><th>Stress (kPa)</th></tr></thead><tbody><tr><td>Scapozza and Bartelt (2003b)</td><td>RG</td><td>[-18.9, -2.1]</td><td>[205, 429]</td><td>1.1 × 10<sup>-6</sup>,<br>4.2 × 10<sup>-6</sup>,<br>1.1 × 10<sup>-5</sup>,<br>4.4 × 10<sup>-5</sup></td><td>[1.7, 258]</td></tr><tr><td>Theile et al. (2011) S1 S2</td><td>RG</td><td>-11</td><td>[205, 345]</td><td>[10<sup>-7</sup>, 10<sup>-6</sup>]</td><td>[2, 10]</td></tr><tr><td>Schleef et al. (2014b)</td><td>PP</td><td>[-18, -3]</td><td>[55, 120]</td><td>[10<sup>-6</sup>, 10<sup>-5</sup>]</td><td>0.133, 0.215, 0.318</td></tr><tr><td>Bernard (2023) V1</td><td>DF/RG</td><td>-18.5 ± 0.5</td><td>[218, 467]</td><td>1.8 × 10<sup>-6</sup></td><td>[0, 300]</td></tr><tr><td>Bernard (2023) V3 V4</td><td>DF/RG</td><td>-18.5 ± 0.5</td><td>[225, 420]</td><td>3.1 × 10<sup>-5</sup></td><td>[0, 300]</td></tr><tr><td>Bernard et al. (2023)</td><td>DF</td><td>-8 ± 0.5</td><td>[236, 260]</td><td>[10<sup>-8</sup>, 10<sup>-7</sup>]</td><td>2.1</td></tr></tbody></table>

The type of snow is defined according to Fierz et al. (2009). The numbers shown in brackets, e.g., [a, b], mean that the parameter ranges between, e.g., a and b. Details can be found in the associated references.

ponent $n(\Phi)$ is the one that minimizes the scatter $\sigma_0(n, \Phi)$ for a given narrow range of $\Phi$ across multiple experiments. We will use this conceptual procedure to derive the stress exponent from the presented experimental data (Sect. 3.2). This problem can be formulated as the minimization of the cost function $\mathcal{L}$ (model parameters) defined as:

$$ \mathcal{L}(\text{model parameters}) = \frac{1}{N_{\text{dataset}}} \sum_{i \in \text{dataset}} \frac{(\tilde{\sigma}_0(\Phi_i) - \sigma_{0i})^2}{\sigma_{0i}^2} \quad (6) $$

and where $\sigma_{0i}$ and $\tilde{\sigma}_0$ denote the measured and model-predicted reference stress, respectively.

## 2.3 Snow settlement in detailed snowpack models

The settlement laws currently implemented in snowpack models such as *Crocus* (Brun et al., 1992; Vionnet et al., 2012) and *Snowpack* (Lehning et al., 2002; Simson et al., 2021) (for vertical stresses lower than 0.4 MPa) are linear. For a given snow layer of thickness $D$, the vertical compaction under an overburden stress $\sigma$ is expressed as:

$$ \frac{dD}{D} = -\frac{\sigma}{\eta} dt \quad (7) $$

where $dr$ is the model time step and $\eta$ (in MPa s) denotes the compaction viscosity. The compaction viscosity is typically parametrised as a function of temperature $T$ (in °C) and density $\rho$ (in kg m$^{-3}$).

In the *Crocus* model (Vionnet et al., 2012), for dry, non-faceted snow, the compaction viscosity follows the formulation of Brun et al. (1992) (hereafter referred to as BR92), derived from the work of Navarre (1975):

$$ \eta_{\text{BR92}} = 4\eta_0 \frac{\rho}{c_\eta} \exp(-a_\eta T + b_\eta \rho) \quad (8) $$

where $\eta_0 = 7.62237$ MPa, $a_\eta = 0.1$ °C$^{-1}$, $b_\eta = 0.023$ m$^3$ kg$^{-1}$, and $c_\eta = 250$ kg m$^{-3}$.

The ensemble version of the French snowpack model *Crocus* (ESCROC) (Lafaysse et al., 2017) also includes

the viscosity parameterisation of Teufelsbauer (2011) (hereafter referred to as T11), derived from the observations of Narita (1984) and corrected for low-density snow using the laboratory-based model of Abe (2001):

$$ \eta_{\text{T11}} = 5 \times 10^{-8} \rho^{(-0.0371T + 4.4)} (10^{-4} \exp(0.018\rho) + 1) \quad (9) $$

In the *Snowpack* model (Lehning et al., 2002), for dry new snow (with dendricity greater than zero), the compaction viscosity follows the empirical formulation of Kojima (1975) (hereafter referred to as K75), derived from field observations of the natural settlement of seasonal snow:

$$ \eta_{\text{K75}} = 7 \times 10^{-9} \rho^{(4.75-T/40)} \quad (10) $$

Although viscosity is a material parameter that characterises deformation in a specific direction, it is commonly assumed that for densities below 450 kg m$^{-3}$, lateral boundary conditions have a negligible effect (Wautier et al., 2017) because of a very small viscous Poisson ratio. Consequently, free and confined boundary conditions are not generally distinguished in these models.

# 3 Results and discussion

## 3.1 Microstructure-based model

Figure 5 presents the results from the microstructure-based simulations conducted on the 37 tested samples.

Figure 5a shows the evolution of the stress exponent $n$ with the solid fraction $\Phi$. The stress exponent ranges between 2.1 and 2.2 for the considered microstructures with $\Phi$ between 0.1 and 0.6. These values are lower than the exponent of dense ice ($n_{\text{ice suquet}} = 2.94$) and are close to the basal slip exponent ($n_{\text{basal}} = 2$). As shown by Védrine et al. (2025) on synthetic polycrystals, the homogenized exponent depends on the activity of the slip systems. Specifically, as the incompatibility of deformation between crystals decreases, the activity of the softer slip systems increases, and the stress
<figure><img src="image_9.png" /><figcaption>Figure 5. Parameters of the homogenized model: (a) stress exponent n and (b) reference stress σ₀ as a function of the solid fraction Φ.</figcaption></figure>

exponent n tends towards 2. We do not observe a large sensitivity of n to Φ, it is rather constant with an average value of nµCT = 2.15. Indeed, the surface ratio between the intercrystalline surfaces and the free surfaces is almost constant and remains small (approximately 0.1, see Védrine and Hagenmuller (2025) and Fig. B1 in Appendix B), which yields a low frustration of the basal slip systems (Védrine et al., 2025). Consequently, we do not observe an increase of n with Φ.

The reference stress σ₀ increases with the solid fraction Φ (Fig. 5b). Indeed, when the solid fraction increases, the stress supported by the ice skeleton (for a given macroscopic stress σ averaged in air and ice) decreases, leading to an increase of the reference stress. Following Védrine et al. (2025), we fitted the evolution of σ₀ with Φ with a power-law:

$$ \sigma_0 = \sigma_0(\Phi = 1)\Phi_r^m. \qquad (11) $$

with Φ_r the rescaled solid fraction (Green, 1998; Holman and Leuenberger, 1988; Stauffer and Aharony, 2018) defined as:

$$ \Phi_r = \frac{\Phi - \Phi_t}{1 - \Phi_t} \qquad (12) $$

The quantity σ₀(Φ = 1) = 272 MPa represents the reference stress of dense polycrystalline ice, m is the solid fraction exponent, and Φ_t is the percolation threshold. By fitting the data obtained through homogenization, we found Φ_t ≈ 0, m = 3.96, and a scatter of L = 1.68 (Eq. 6) around this model. The very low value for Φ_t suggests that the snow microstructure can maintain connectivity (i.e., percolate) even at very low solid fractions. The solid fraction sensitivity m relates to the heterogeneity of viscoplastic deformations within the sample. More generally, the value of m can be computed for different snow properties, such as thermal conductivity,

Young's modulus, and strength. It serves as an indicator of the heterogeneity of the local driving processes: higher values of m reflect higher heterogeneity (modulated by both the process type and microstructure) (Bruno et al., 2011). The fitted value for snow viscoplasticity is of m = 3.96. For snow thermal conductivity, m typically ranges between 1 and 2 (Calonne et al., 2011, 2019). For the elastic modulus, m usually falls between 2.5 and 3.5 (Gaume et al., 2017; Sundu et al., 2024a; Blatny et al., 2021). For strength, Ritter et al. (2020) found values of 3.59 in compression and 3.09 in tension. To our knowledge, a universal relationship between microstructure morphology and the m value for a given physical property is still unknown.

Overall, the pseudo-linear viscosity η scales with φ⁻¹²m (Eq. 2). This extreme sensitivity to solid fraction emphasizes the challenges in characterizing the viscoplastic behaviour of snow. For example, a 5% relative error in measuring the solid fraction, which is common (e.g. Proksch et al., 2015; Hagenmuller et al., 2016) leads to a relative error of about 20% error in the reference stress σ₀ and a relative error of around 50% for the compaction viscosity η.

The remaining scatter (L = 1.68) can be attributed to the fact that only density is considered as a predictive indicator (Fig. 5b). For instance, RGlr exhibits a reference stress approximately four times lower than that of other snow types at the same density. Similarly, the reference stress simulated for the DF series studied by Bernard (2023) cannot be perfectly fitted with a constant exponent m. We observe that the values of m tend to decrease with the densification. These two examples suggest the influence of additional geometric parameters, such as bond size (Hagenmuller et al., 2014a), which are currently not accounted for in the model.
## 3.2 Model calibrated on mechanical tests

To reinterpret previous experimental data, our idea was to determine the value of *n* that minimizes the scatter of the reference stress σ₀ for “similar” solid fractions Φ across different experiments (see Sect. 3.1 and Eq. 11). The microstructure-based simulations showed that the reference stress is likely to be very sensitive to Φ: σ₀ ~ Φ³.⁹⁶. It is thus impossible from the collected data (178 points) to gather enough points within a narrow range of [Φ, Φ+ΔΦ] with ΔΦ ≪ 1 to accurately quantify the scatter of σ₀ and the value *n* (Φ). However, the microstructure-based simulations also revealed that *n* is relatively constant for snow and σ₀ follows a power-law of Φ. We use these two results from the microstructure-based simulations as additional assumptions to reinterpret previous experimental data. In other words, we want to find a constant *n* that minimizes the relative scatter *L(n)* (defined as the relative standard deviation, see Eq. 6) of σ₀ around a power-law of Φ of the form $\tilde{\sigma}_0(n) = \sigma_0(\Phi = 1)\Phi_r^m$, with $\Phi_r$ the rescaled solid fraction (Eq. 12) and, *m* and $\Phi_t$ two optimal fitting parameters depending on *n*.

We tested this procedure with three different values for the stress exponent: *n* = 1, *n* = 3, and *n* = 2.15 (Fig. 6):

- *n* = 1 (Fig. 6a). When the viscoplastic behaviour is assumed to be linear, we observe significant dispersion in the reference stress σ₀ (*L* = 0.87), particularly at solid fractions exceeding 0.2. Notably, the dataset is divided into two distinct groups: the tests with imposed strain rates are located in the lower part of the scatter plot, while the creep tests are located in the upper part. This dispersion and the separation into two groups indicate that a linear model fails to adequately reproduce the variety of experimental data. Figure 6 also includes the linear models from Brun et al. (1992), Teufelsbauer (2011) and Kojima (1975) at a temperature of -10 °C. The Br92 model overestimates all experimental reference stresses and thus appears too fluid. Moreover, it reaches the viscosity of ice for an ice volume fraction of about 0.5. The T11 model also tends to overestimate the reference stress but aligns more closely with the experimental values for low-density snow measured by Schleef et al. (2014b). Finally, the K75 model reproduces the experiments of Schleef et al. (2014b) with good fidelity and lies near the centre of the experimental data set.

- *n* = 3 (Fig. 6b). The case *n* = 3 corresponds to the exponent of dense polycrystalline ice. It is commonly assumed that this exponent of the dense material applies to its porous counterpart, snow (Kirchner et al., 2001; Wautier et al., 2017; Fourteau et al., 2023). In this case, the dispersion decreases to *L* = 0.19. This decrease is particularly visible for snow with solid fractions above 0.2. However, the values of σ₀ are still divided into two groups depending on the type of test.

- *n* = 2.15 (Fig. 6c). Finally, when we assume that the stress exponent is equal to *n*<sub>μCT</sub> = 2.15, the dispersion is reduced to *L* = 0.14, with optimized parameters Φ<sub>t</sub> = 0.025 and *m* = 3.195. Moreover, unlike the case with *n* = 3, the experiments with imposed stress and strain rates are superposed, which indicates an appropriate underlying model.

Figure C1 shows the evolution of the dispersion *L* as a function of *n*. *L* decreases when *n* increases from 1 to 2.35. For *n* in [2.35, 2.59], *L* remains nearly constant with a weak minimum at *n* = 2.53. For *n* > 2.59, *L* increases with *n*. Interestingly, the minimum dispersion error is obtained for a value of *n* close to the one obtained from microstructure-based simulations. This parametrization significantly improves the modelling of the viscoplastic behaviour of snow. The exponent *n* = 2.15 is consistent with the value derived by Scapozza and Bartelt (2003b) and Schleef et al. (2014b) on samples with low density (Φ < 0.2). This agreement supports the idea that basal glide is the dominant mechanism at low density. Furthermore, this suggests that the reduction of the stress exponent from dense ice to snow does not require introducing grain boundary sliding (GBS) as an additional deformation mechanism. GBS was originally proposed by Alley (1987) and Goldsby and Kohlstedt (1997) to explain observations of flow stress sensitivity at low solid fractions and is associated with a stress exponent of *n* ≈ 2 in fine-grained ice (Goldsby and Kohlstedt, 2001). Our work shows that the low crystalline frustration in snow provides an alternative physical explanation for the reduced stress exponent observed in snow.

The residual error obtained, *L(n* = 2.15) = 0.15 is not null and may be attributed to model error and experimental measurement noise. First, the microstructure is simplified here to its solid fraction, Φ. Our model cannot account for other geometrical properties, such as bond size (Hagenmuller et al., 2014a) or connectivity (Schleef et al., 2014b; Schöttner et al., 2025), or the inter-crystalline surface area, which affects the frustration of the basal glide (Védrine et al., 2025). Indeed, the residual error computed on the microstructure-based simulations is also not null (*L* = 1.68), which tends to show that the variety of microstructural patterns for a given Φ could be enriched by other microstructural descriptors. Last, we showed that σ₀ varies like Φ<sup>m</sup> with *m* ≃ 3. Therefore, a small error in measuring Φ, or a small heterogeneity of density in the prepared sample, may significantly affect the measured points and subsequent analysis. Second, the scaling to a constant temperature of -10 °C might be a source of uncertainty. For instance, the residual dispersion between appears more pronounced at solid fractions of 0.4 and 0.5 (Fig. 6c), corresponding to tests conducted by the same operator (Scapozza and Bartelt, 2003a) under the same conditions except different temperatures ([~ -18.5, ~ -11, ~ -2]). Likely, the Arrhenius scaling around -10 °C (Eq. 3) does not apply to temperatures close to 0 °C and below -15 °C, where different
<figure><img src="image_10.png" /><figcaption>**Figure 6.** Experimental data on the evolution of the reference stress as a function of solid fraction, recalculated under various parametrisations of *n* at a temperature of -10 °C.</figcaption></figure>

deformation mechanisms may take place (Mellor and Testa, 1969b).

Figure 7 illustrates the reference stress $\sigma_0$ computed from microstructure-based simulations alongside the experimental re-analysis, both with $n = 2.15$. We observe that the microstructure-based reference stress deviates from the experimental reference stress by a factor of 2. This difference indicates that the viscoplastic behaviour simulated by

considering only intracrystalline deformations might be too fluid. Thus, incorporating an additional deformation mechanism, such as grain boundary sliding, would likely render the behaviour even more fluid. One possible explanation for this lower viscosity may be the assumption that a geometric grain corresponds to a crystal (Védrine et al., 2025). Finally, both fits to experimental data and microstructure-based simulations tend to underestimate viscosity for a solid fraction
<figure><img src="image_11.png" /><figcaption>**Figure 7.** Comparison of the microstructure-based model and experimental data regarding the evolution of the reference stress as a function of solid fraction, calculated with the assumption $n = 2.15$ at a temperature of $-10$ °C.</figcaption></figure>

around 0.4, suggesting a potential evolution of the morphology parameter $m$, which is here assumed constant for $\Phi$ in the range [$\Phi_t$, 1]. Indeed, Coble (1961) showed that for $\Phi$ close to 1, i.e. diluted pores in an ice matrix, the solid fraction exponent $m$ is closer to 1.

### 3.3 Difference between the non-linear model calibrated on mechanical tests and the settlement law implemented in Crocus

In this section, we compare the predictions of the detailed snowpack model Crocus (referred to as Br92) to those of our non-linear model calibrated on mechanical tests (referred to as V25):

$$ \dot{\epsilon}(\sigma, T) = e^{-\frac{Q}{R}\left(\frac{1}{T} - \frac{1}{T_0}\right)} \left(\frac{\sigma}{\sigma_0(\Phi=1)\Phi_r^m}\right)^n, \quad (13) $$

with $Q = 69.1 \text{ kJ mol}^{-1}$, $T_0 = 263$ K, $n = 2.15$, $m = 3.195$ and $\Phi_r$ the rescaled solid fraction (Eq. 12) with $\Phi_t = 0.025$. The evaluation is focused on stress - solid fraction ($\sigma, \Phi$) couples, which are realistic for natural snowpacks. We recall that the temporal evolution of the thickness $h$ of a layer subjected to an overburden stress $\sigma$ is related to strain rate as $\dot{\epsilon} = \dot{h}/h$. Note that the overburden stress is defined as the gravitational acceleration ($g = 9.81 \text{ m s}^{-2}$) multiplied by the snow water equivalent of the overlying snow.

Figure 8 shows the relationship between the applied stress and the solid fraction for layers measured in natural seasonal snowpacks at different locations: Cambridge Bay (High-Arctic, Canada, winter 2023–2024, ongoing ERC IVORI campaign), Col de Porte (French Alps, 1325 m, winters 1993–1994 to 2017–2018, Lejeune et al., 2019), and Weissflüjoch (Swiss Alps, 2536 m, winter 2015–2016, Calonne et al., 2020). The applied stress was calculated as the overburden at mid-layer. It mainly ranges between 0.1 and 4 kPa and generally increases with the solid fraction. The two alpine sites exhibit larger stresses for a given solid fraction compared to the Arctic site. In the dry arctic region, densification is primarily driven by wind erosion-deposition, whereas at the alpine sites, gravity-driven creep is the dominant settlement mechanism (Woolley et al., 2024).

The predictions of V25 and Br92 are compared in the stress – solid fraction diagram shown in Fig. 8. The Br92 model exhibits a lower stress sensitivity due to its linear formulation. Consequently, our parametrisation V25 predicts higher strain rates for large loads and lower strain rates for small loads, compared to Br92. Notably, the two models yield equivalent results for the solid fraction-stress couples most commonly observed at the Col de Porte, where the Br92 model has been extensively calibrated (Navarre, 1975; Brun et al., 1992; Vionnet et al., 2012; Lafaysse et al., 2017). This agreement indicates that V25 behaves similarly to Br92 under the natural conditions encountered on this specific site, and therefore suggests that the V25 model is also likely to reproduce the snow compaction observed at Col de Porte. The linear model performs well when the applied stress does not vary much beyond its calibration points or when the stress correlates with density. In this latter case, the non-linearity can be artificially introduced through the dependence of viscosity $\eta$ on solid fraction $\Phi$. However, substantial discrepancies arise when moving beyond these “standard” conditions: Br92 predicts higher strain rates than V25 for $\Phi > 0.08$ and stresses smaller than those observed at Col de Porte. Conversely, for very low density samples ($\Phi < 0.08$) and for stresses higher than the “standard” ones, Br92 predicts a lower strain rate compared to V25.

We illustrated this discrepancy on three exemplary cases, which are out of the common conditions encountered at Col de Porte:

- Heavy snowfalls. We consider the bottom of a 25 cm thick layer of snow with a density of $70 \text{ kg m}^{-3}$ ($\Phi = 0.076$, $\sigma = 0.17 \text{ kPa}$, see green point on Fig. 8). This scenario typically corresponds to a layer formed during a heavy snowfall. During the first day, the V25 model shows densification at a rate up to 33 times that predicted by the Br92 model (Fig. 9). Interestingly, the two models converge after approximately 14 d, as the increase in density leads to a corresponding rise in viscosity. This finding aligns with field observations, which often indicate that densification is frequently underesti-
<figure><img src="image_12.png" /><figcaption>**Figure 8.** Normalized residual $R = \frac{\dot{\epsilon}_{V25} - \dot{\epsilon}_{Br92}}{\epsilon_{Br92} + \epsilon_{V25}}$ between the V25 and Br92 model for varying solid fractions and stresses. The stress is converted into its snow water equivalent.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 9.** Evolution of the density at the bottom of a 25 cm layer of new snow at 70 kg m⁻³ with the model Br92 et V25 at -10 °C.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 10.** Evolution of the density at the bottom of a 10 cm thick layer of snow at 250 kg m⁻³ with the model Br92 et V25 at -10 °C.</figcaption></figure>

mated in the immediate aftermath of snowfall (Lundy et al., 2001; Steinkogler et al., 2009; Wever et al., 2015) although it has minimal impact in the long term. Quéno et al. (2016) reached similar conclusions with AROME–Crocus (included BR92 model) when comparing simulated and observed snow depth. Higher strain rates thus

naturally arise on recent snow in the V25 model without introducing a distinct role for curvature-driven metamorphism (Chen et al., 2019).

- *Wind-packed surface layers.* We consider the bottom of a 10 cm thick layer of snow with a density of 250 kg m⁻³ (Φ = 0.22, σ = 0.25 kPa, see blue point on
<figure><img src="image_15.png" /><figcaption>**Figure 11.** Evolution of the density of a layer of firn at 450 kg m⁻³ buried below 2 m of snow with an average density of 400 kg m⁻³) with the model Br92 et V25 at -10 °C.</figcaption></figure>

Fig. 8). This layer typically forms at the snowpack surface during a wind event. The V25 model predicts a densification rate about six times slower in the first few days compared to the Br92 model (see Fig. 10). This result aligns with observations of limited densification of basal layers in Arctic snowpacks—behaviour that is not well captured by existing parametrisations (Woolley et al., 2024). Nevertheless, our model does not distinguish snow types and cannot predict the distinct densification rates observed between RG and DH at a given density (Fourteau et al., 2023).

- *Deep firn.* We consider a layer of firn with a density of 450 kg m⁻³ buried beneath 2 m of snow with an average density of 400 kg m⁻³ (Φ = 0.49, σ = 11.9 kPa, see yellow point on Fig. 8). This layer typically forms on an accumulation zone on an alpine glacier. The V25 model predicts a densification rate that can be up to 30 times faster than that predicted by the Br92 model (Fig. 9). This result is consistent with the fact that the BR92 model predicts the viscosity of ice at a solid fraction of 0.5. In addition, flow rates of dense snow and firn in glacier accumulation zones are frequently underestimated by Br92 (Touzeau et al., 2018; Verjans et al., 2019). To trick this discrepancy, van Kampenhout et al. (2017) tuned the viscosity of BR92 by a factor of 1.4 to model polar firn settlement. Nevertheless, to broaden the applicability of the V25 model to firn at densities exceeding 450 kg m⁻³, we will need to investigate the effects of lateral confinement in greater detail, which are negligible only for low-density snow (Wautier et al., 2017).

## 3.4 Intrinsic limitations of the proposed settlement model

We have integrated previously separate approaches into a cohesive, physically grounded framework, allowing us to reconcile results that may initially appear contradictory. However, our framework (Eq. 13) cannot, at this stage, capture known features of snow settlement because it over-simplifies the snow microstructure or the underlying processes. The limitations that would need to be addressed in future work notably comprise:

- *Accounting for morphology descriptors.* Our model simplifies the snow microstructure to its ice fraction Φ. This common approach does not capture second-order effects of the microstructure on its mechanical behaviour (Hagenmuller, 2014). Indeed, we observed scatter of the simulated reference stress at a given density (Fig. 5). We chose not to introduce any additional microstructural descriptors to maintain compatibility with existing experimental data or detailed snowpack models in which the microstructure is very roughly described. However, conducting the microstructure-based simulations on a larger set of 3D images could help identify key additional descriptors of the microstructure (e.g. Hagenmuller et al., 2014a; Schöttner et al., 2025). These descriptors will need to be predicted by the next generation of snowpack models (e.g. Brondex et al., 2023). Therefore, monitoring these indicators during in-situ compression tests (e.g. Bernard, 2023) could help identify their evolution laws.

- *Accounting for the inter-crystalline surface area.* Our microstructure-based simulations assume that one grain detected with geometrical criteria (bond ~ neck) corresponds to one crystal. This assumption is reasonable for rounded grains, but is more questionable for faceted crystals or depth hoar. Védrine et al. (2025) showed that the value of the stress exponent *n* is sensitive to crystal frustration, which depends on the ratio between the free surface area and the inter-crystalline surface area. Based on our assumption, the ratio of inter-crystalline to free surface area was found to be nearly constant and small (approximately 0.1), implying a low degree of basal slip-system frustration. We believe that this inter-crystalline surface area plays a key role in the observed differences in settlement rates between rounded grains and depth hoar (Kojima, 1967; Armstrong, 1980; Sturm and Holmgren, 1998). In particular, we intuit that, for a given density, depth hoar is characterised by a higher frustration compared to rounded grains, which yields slower densification. This lower settlement rate of depth hoar is not anecdotal. It is crucial in the formation of persistent weak layers, a prerequisite for slab avalanches, or in capturing the thermal properties of widespread tundra snowpacks (Fourteau et al., 2023).
Due to the lack of knowledge about the exact crystalline structure (i.e., crystal orientation and intercristalline surface) of depth hoar, we cannot currently test this intuition. This idea thus warrants further numerical investigations using the proposed 3D model associated with explicit crystal measurements via diffraction contrast tomography (Granger et al., 2021).

- Refining the temperature sensitivity (Arrhenius law). We assumed that the temperature dependence of the homogenised power's law follows an Arrhenius relationship with a single activation energy applied uniformly across all microstructures and temperatures (Eq. 3). However, the activation energy in dense polycrystalline ice varies with temperature (Mellor and Testa, 1969b), and each crystallographic slip system has its own activation energy (Mellor and Testa, 1969b; Ramseier, 1972; Goldsby and Kohlstedt, 2001), affecting the relative activity of different systems. Consequently, a temperature dependence of the stress exponent is to be expected, as noted by Scapozza and Bartelt (2003b). Besides, explicit measurements of the activation energy for snow are seldom reported (Delmas, 2013; Schleef et al., 2014b) and can be affected by microstructural changes between tests or during a single test. Knowledge of the temperature sensitivity of the driving deformation mechanisms would help develop a unified settlement law applicable over a broad range of temperatures.

- Extending the approach to a 3D constitutive model. We evaluated our model on snow with a solid fraction, $\Phi \le 0.5$, but it could be extended to larger solid fractions, e.g., to capture the transition from snow to glacier ice. However, our model is one-dimensional and assumes uniaxial compression with zero average horizontal stresses. Snow, firn or ice in natural conditions may also be subjected to horizontal stresses. For snow with a low solid fraction, these stresses can generally be neglected (Wautier et al., 2017), but this is not the case any more for firn or porous ice. Therefore, it would be necessary to develop a complete 3D constitutive model to capture this effect.

# 4 Conclusions and perspectives

We used three different sources of information to understand the viscoplastic behaviour of snow:

- Microstructure-based simulations. We simulated the viscoplastic behaviour of snow based on its three-dimensional microstructure and a crystal-plasticity law for the ice matrix. Our simulations revealed, across a variety of microstructures, that the viscoplastic behaviour of snow followed a power-law $\dot{\epsilon} = \dot{\epsilon}_0(\sigma/\sigma_0)^n$. The stress exponent $n$ ranged from 2.1 to 2.2. The

reference stress $\sigma_0$ was found to be a function of the solid fraction $\Phi$, following the relationship $\sigma_0 \sim \Phi_m^n$. The values for $n$ are consistent both with experimental results for low-density snow (Scapozza and Bartelt, 2003b; Schleef et al., 2014b) and with generic predictions for porous polycrystals (Védrine et al., 2025). These values are close to the stress exponent of the basal dislocation glide, since the ice crystals are poorly constrained by neighbours in snow.

- Previous viscoplastic tests. We revisited a dataset comprising 178 measurement points $(\sigma, \dot{\epsilon})$ and found that linear models fail to reproduce the viscoplastic behaviour of snow under various loading conditions. In contrast, using the stress exponent derived from microstructure-based simulations ($n = 2.15$) significantly reduces the variability of the reference stress between independent measurement sets, especially between load-controlled and displacement-controlled tests. Additionally, it was unnecessary to adjust $n$ over a wide range (e.g., [1, 4], which exceeds the range predicted by the microstructure-based model [2, 3]) to consolidate previous datasets. We obtained a relatively simple constitutive relation presumably valid across a broad range of loading conditions (Eq. 13).

- Detailed snowpack models. Snowpack models, such as Crocus, assume a linear relation between stress and strain rate ($n = 1$) and are well-calibrated on field data. We showed that in typical natural seasonal alpine snowpacks, the applied stress (or layer depth) is not independent of the ice fraction. Therefore, linear and non-linear models can coincide under these specific conditions. In fact, the linear viscosity $\eta$, in snowpack models, implicitly accounts for both the material viscosity variation with $\Phi$ and the site-specific correlation between the stress and $\Phi$. As a result, $\eta$ can artificially introduce non-linearity into the constitutive equation. Indeed, we showed that Crocus and our model yield equivalent predictions on conditions typically encountered on alpine sites. However, large deviations appear under “non-standard” conditions, such as high loads on light snow or low loads on dense snow, where our model appears to perform better qualitatively.

We have thus unified previously disparate approaches within a single, physically grounded framework, thereby reconciling results that might initially seem contradictory. This framework provides a robust basis for improving the representation of snow settlement processes in numerical snowpack models.

# Appendix A: Constitutive Law: Crystal plasticity model

The total strain rate tensor $\epsilon$ is decomposed into an elastic part ($\epsilon^e$) and a viscoplastic part ($\epsilon^{vp}$), such that: $\epsilon = \epsilon^e + \epsilon^{vp}$.
<table><caption>**Table A1.** Single crystal ice parameters of crystal plasticity models at -10 °C (Suquet et al., 2012).</caption><thead><tr><th>Family</th><th>Systems</th><th>n<sup>(k)</sup></th><th>τ<sub>0</sub><sup>(k)</sup> (MPa)</th></tr></thead><tbody><tr><td>Basal [0001] &lt; 11&#x0304;20 &gt;</td><td>3</td><td>2</td><td>35</td></tr><tr><td>Prismatic [01&#x0304;10] &lt; 2&#x0304;10 &gt;</td><td>3</td><td>2.85</td><td>166</td></tr><tr><td>Pyramidal [11&#x0304;22] &lt; 11&#x0304;23 &gt;</td><td>6</td><td>4</td><td>98</td></tr></tbody></table>

The elastic strain $\varepsilon^e$ is related to the stress tensor $\sigma$ and the fourth-order stiffness tensor $C$ by the relation:

$$ \sigma = C : \varepsilon^e. \qquad (A1) $$

The elasticity tensor $C$ is assumed to be isotropic transverse in the crystal reference frame and is defined according to Gammon et al. (1983) as follows: $C_{11} = 13.9$ GPa, $C_{33} = 15.0$ GPa, $C_{44} = 3.0$ GPa, $C_{12} = 7.1$ GPa, $C_{13} = 5.8$ GPa.

For the viscoplastic part, at infinitesimal strains, the vis-
coplastic strain results from slips on a total of $N$ different
slip systems:

$$ \dot{\varepsilon}^{vp} = \sum_{k=1}^{N} \dot{\gamma}_{0}^{(k)}(T) \left( \frac{|\tau^{(k)}|}{\tau_{0}^{(k)}} \right)^{n^{(k)}} \operatorname{sgn}(\tau^{(k)}) \mu^{(k)} \quad (A2) $$

where $n^{(k)}$ is the stress exponent, $\dot{\gamma}_0^{(k)}$ is the reference shear rate, $\tau_0^{(k)}$ is the critical shear stress for slip system (k), and $\mu^{(k)} = n^{(k)} \otimes_s b^{(k)}$ is the Schmid tensor defined by the (symmetric) dyadic product $(\otimes_s)$ of the vector normal to the glide plane $n^{(k)}$ and the Burgers vector $b^{(k)}$ for slip system (k). As we only focus on the steady-state response and to simplify comparison, the critical stresses $\tau_0^{(k)}$ are assumed to be constant with time (no hardening), and the reference shear rate is set to $\dot{\gamma}_0^{(k)} = 1 \text{ s}^{-1}$. We calibrate the values of the Critical Resolved Shear Stresses (CRSS) to obtain the viscoplastic law of dense polycrystalline ice proposed by Budd and Jacka (1989); Castelnau et al. (1996) at -10 °C ($A = 7.8 \times 10^{-8}$ MPa$^{-1}$ s$^{-3}$ and $n_{\text{poly}} = 2.94$ the parameters of the Glen law for polycrystalline ice $\dot{\epsilon} = A\sigma^{n_{\text{poly}}}$). The parametrizations of the crystal plasticity model are summarised in Table A1.

## Appendix B: Evolution of the surface ratio between the inter-crystalline surfaces and the free surfaces with the microstructure

The evolution of the surface ratio between inter-crystalline surfaces and free surfaces, denoted as $r$, is plotted as a function of the solid fraction for different microstructures (Fig. B1). We observe that $r$ increases only slightly with the solid fraction. A plateau is reached at $r = 0.1$ for solid fractions ranging from 0.2 to 0.45, with only the MM and RGlr snow types deviating from this trend. This low surface ratio leads to low frustration of the basal slip systems (Védrine et al., 2025). Note that the surface ratio obtained for these snow microstructures is significantly lower than that for Voronoi microstructures with pores and crystals of similar size (Védrine et al., 2025). Indeed, for a solid fraction of $\phi = 0.4$, the snow microstructures exhibit a surface ratio of approximately $\beta \approx 0.1$, whereas the Voronoi microstructures exhibit a value of about $\beta \approx 0.35$.

<figure><img src="image_16.png" /><figcaption>**Figure B1.** Contiguity $\beta$ and surface ratio $r$ as the function of solid fraction $\Phi$. The contiguity Underwood (1970) is defined as $\beta = \frac{r}{1+r}$.</figcaption></figure>
**Appendix C: Evolution of $\mathcal{L}$ as a function of stress exponent $n$**.

<figure><img src="image_17.png" /><figcaption>**Figure C1.** Evolution of the dispersion $\mathcal{L}$ as a function of stress exponent $n$ at a temperature of -10 °C.</figcaption></figure>

*Code and data availability.* All materials used in this article (codes, geometric characteristic of the microstructure, numerical values, experimental data etc.) are available at https://doi.org/10.5281/zenodo.16778695 (Védrine and Hagenmuller, 2025).

*Author contributions.* LV developed the microstructure-based model, performed the simulations, compiled the data, and created the figures. LV and PH interpreted the results and wrote the paper. PH designed the study.

*Competing interests.* The contact author has declared that neither of the authors has any competing interests.

*Disclaimer.* Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. The authors bear the ultimate responsibility for providing appropriate place names. Views expressed in the text are those of the authors and do not necessarily reflect the views of the publisher.

*Acknowledgements.* The computations presented in this paper were performed using the GRICAD infrastructure (https://gricad.univ-grenoble-alpes.fr, last access: 13 February 2026), which is supported by CNRS, University Grenoble Alpes and INRIA. We thank J. Schweizer, and O. Castelnau and two anonymous reviewers for their constructive feedback on this work. We thank H. Löwe for fruitful discussions.

*Financial support.* CEN is part of LabEx OSUG@2020 (ANR-15-IDEX-02). We thank the European Research Council (ERC) for funding, under the European Union's Horizon 2020 research and innovation programme (grant agreement no. 949516, IVORI), the Arctic field campaign in Cambridge Bay. The tomography apparatus (TomoCold) was funded by INSU-LEFE, Labex OSUG (Investissements d'avenir-ANR10 LABX56).

*Review statement.* This paper was edited by Jürg Schweizer and reviewed by Olivier Castelnau and one anonymous referee.

**References**

Abe, O.: Creep experiments and numerical simulations of very light artificial snowpacks, Annals of Glaciology, 32, 39–43, https://doi.org/10.3189/172756401781819201, 2001.

Alley, R. B.: Firn densification by grain-boundary slidinf: a first model, Le Journal de Physique Colloques, 48, C1–C256, https://doi.org/10.1051/jphyscol:1987135, 1987.

Anderson, D. G.: Iterative Procedures for Nonlinear Integral Equations, Journal of the ACM, 12, 547–560, https://doi.org/10.1145/321296.321305, 1965.

Armstrong, R. L.: An analysis of compressive strain in adjacent temperature-gradient and equi-temperature layers in a natural snow cover, Journal of Glaciology, 26, 283–289, https://doi.org/10.3189/S0022143000010820, 1980.

Arnaud, L., Lipenkov, V., Barnola, J. M., Gay, M., and Duval, P.: Modelling of the densification of polar firn: characterization of the snow-firn transition, Annals of Glaciology, 26, 39–44, https://doi.org/10.3189/1998AoG26-1-39-44, 1998.

Auriault, J.-L., Boutin, C., and Geindreau, C.: Homogenization of Coupled Phenomena in Heterogenous Media, John Wiley & Sons, ISBN 978-0-470-61044-2, 2010.

Bader, H.: Sorge's Law of Densification of Snow on High Polar Glaciers, Journal of Glaciology, 2, 319–323, https://doi.org/10.3189/S0022143000025144, 1954.

Bader, H.: Theory of densification of dry snow on high polar glaciers, II, Report, Cold Regions Research and Engineering Laboratory (U.S.), 1962.

Bahaloo, H., Forsberg, F., Lycksam, H., Casselgren, J., and Sjödahl, M.: Material mapping strategy to identify the density-dependent properties of dry natural snow, Applied Physics A, 130, 141, https://doi.org/10.1007/s00339-024-07288-y, 2024.

Barnett, T. P., Pierce, D. W., Hidalgo, H. G., Bonfils, C., Santer, B. D., Das, T., Bala, G., Wood, A. W., Nozawa, T., Mirin, A. A., Cayan, D. R., and Dettinger, M. D.: Human-Induced Changes in the Hydrology of the Western United States, Science, 319, 1080–1083, https://doi.org/10.1126/science.1152538, 2008.

Barnola, J. M., Raynaud, D., Korotkevich, Y. S., and Lorius, C.: Vostok ice core provides 160,000-year record of atmospheric CO2, Nature, 329, 408–414, https://doi.org/10.1038/329408a0, 1987.

Barraclough, T. W., Blackford, J. R., Liebenstein, S., Sandfeld, S., Stratford, T. J., Weinländer, G., and Zaiser, M.: Propagating compaction bands in confined compression of snow, Nature Physics, 13, 272–275, https://doi.org/10.1038/nphys3966, 2017.
Bartelt, P. and Christen, M.: A computational procedure for instationary temperature-dependent snow creep, in: Advances in Cold-Region Thermal Engineering and Sciences, edited by: Hutter, K., Wang, Y., and Beer, H., 367–386, Springer, Berlin, Heidelberg, ISBN 978-3-540-48410-3, https://doi.org/10.1007/BFb0104195, 1999.

Bergen, J. D.: Some measurments of air permeability in mountain snow cover, International Association of Scientific Hydrology Bulletin, 13, 5–13, https://doi.org/10.1080/02626666809493602, 1968.

Bernard, A.: Etude multiéchelle de la transition ductile-fragile dans la neige, phdthesis, Université Grenoble Alpes, PhD thesis, 2023GRALI027, https://theses.hal.science/tel-04145610v1/file/BERNARD_2023_archivage.pdf (last access: 13 February 2026), 2023.

Bernard, A., Hagenmuller, P., Montagnat, M., and Chambon, G.: Disentangling creep and isothermal metamorphism during snow settlement with X-ray tomography, Journal of Glaciology, 69, 899–910, https://doi.org/10.1017/jog.2022.109, 2023.

Blatny, L., Löwe, H., Wang, S., and Gaume, J.: Computational micromechanics of porous brittle solids, Computers and Geotechnics, 140, 104284, https://doi.org/10.1016/j.compgeo.2021.104284, 2021.

Brondex, J., Fourteau, K., Dumont, M., Hagenmuller, P., Calonne, N., Tuzet, F., and Löwe, H.: A finite-element framework to explore the numerical solution of the coupled problem of heat conduction, water vapor diffusion, and settlement in dry snow (IvoriFEM v0.1.0), Geoscientific Model Development, 16, 7075–7106, https://doi.org/10.5194/gmd-16-7075-2023, 2023.

Brun, E., David, P., Sudul, M., and Brunot, G.: A numerical model to simulate snow-cover stratigraphy for operational avalanche forecasting, Journal of Glaciology, 38, 13–22, https://doi.org/10.3189/S0022143000009552, 1992.

Bruno, G., Efremov, A. M., Levandovskyi, A. N., and Clausen, B.: Connecting the macro- and microstrain responses in technical porous ceramics: modeling and experimental validations, Journal of Materials Science, 46, 161–173, https://doi.org/10.1007/s10853-010-4899-0, 2011.

Budd, W. F. and Jacka, T. H.: A review of ice rheology for ice sheet modelling, Cold Regions Science and Technology, 16, 107–144, https://doi.org/10.1016/0165-232X(89)90014-1, 1989.

Calonne, N., Flin, F., Morin, S., Lesaffre, B., du Roscoat, S. R., and Geindreau, C.: Numerical and experimental investigations of the effective thermal conductivity of snow, Geophysical Research Letters, 38, https://doi.org/10.1029/2011GL049234, 2011.

Calonne, N., Milliancourt, L., Burr, A., Philip, A., Martin, C. L., Flin, F., and Geindreau, C.: Thermal Conductivity of Snow, Firn, and Porous Ice From 3-D Image-Based Computations, Geophysical Research Letters, 46, 13079–13089, https://doi.org/10.1029/2019GL085228, 2019.

Calonne, N., Richter, B., Löwe, H., Cetti, C., ter Schure, J., Van Herwijnen, A., Fierz, C., Jaggi, M., and Schneebeli, M.: The RHOSSA campaign: multi-resolution monitoring of the seasonal evolution of the structure and mechanical stability of an alpine snowpack, The Cryosphere, 14, 1829–1848, https://doi.org/10.5194/tc-14-1829-2020, 2020.

Castelnau, O., Duval, P., Lebensohn, R. A., and Canova, G. R.: Viscoplastic modeling of texture development in polycrystalline ice with a self-consistent approach: Comparison with bound

estimates, Journal of Geophysical Research: Solid Earth, 101, 13851–13868, https://doi.org/10.1029/96JB00412, 1996.

Castelnau, O., Duval, P., Montagnat, M., and Brenner, R.: Elastoviscoplastic micromechanical modeling of the transient creep of ice, Journal of Geophysical Research: Solid Earth, 113, https://doi.org/10.1029/2008JB005751, 2008.

Chandel, C., Srivastava, P. K., and Upadhyay, A.: Estimation of Rheological Properties of Snow Subjected to Creep, Defence Science Journal, 57, 393–401, https://doi.org/10.14429/dsj.57.1786, 2007.

Chen, Y., Gélébart, L., Chateau, C., Bornert, M., Sauder, C., and King, A.: Analysis of the damage initiation in a SiC/SiC composite tube from a direct comparison between large-scale numerical simulation and synchrotron X-ray micro-computed tomography, International Journal of Solids and Structures, 161, 111–126, https://doi.org/10.1016/j.ijsolstr.2018.11.009, 2019.

Chevy, J.: Viscoplasticité et Hétérogénéités de déformation du monocristal de glace: expériences et simulations, phdthesis, Institut National Polytechnique de Grenoble – INPG, https://theses.hal.science/tel-00396410v1/file/2008-these-Juliette-CHEVY.pdf (last access: 13 February 2026), 2008.

Coble, R. L.: Sintering Crystalline Solids, I. Intermediate and Final State Diffusion Models, Journal of Applied Physics, 32, 787–792, https://doi.org/10.1063/1.1736107, 1961.

Colbeck, S. C.: Air Movement in Snow Due to Wind-pumping, Journal of Glaciology, 35, 209–213, https://doi.org/10.3189/S0022143000004524, 1989.

Delmas, L.: Influence of snow type and temperature on snow viscosity, Journal of Glaciology, 59, 87–92, https://doi.org/10.3189/2013JoG11J231, 2013.

Desrues, J., Darve, F., Flavigny, E., Navarre, J. P., and Taillefer, A.: An Incremental Formulation of Constitutive Equations for Deposited Snow, Journal of Glaciology, 25, 289–307, https://doi.org/10.3189/S0022143000010509, 1980.

Dormieux, L. and Bourgeois, E.: Introduction to porous media micro-mechanics; Introduction a la micromecanique des milieux poreux, Presses de l'Ecole nationale des ponts et chaussées, https://www.osti.gov/etdeweb/biblio/20479238 (last access: 13 February 2026), 2002.

Duval, P., Ashby, M. F., and Anderman, I.: Rate-controlling processes in the creep of polycrystalline ice, The Journal of Physical Chemistry, 87, 4066–4074, https://doi.org/10.1021/j100244a014, 1983.

Fierz, C., Armstrong, R. L., Durand, Y., Etchevers, P., Greene, E., McClung, D. M., Nishimura, K., Satyawali, P. K., and Sokratoğlu, S. A.: The international classification for seasonal snow on the ground, IHP-VII Technical Documents in Hydrology no. 83, IACS Contribution no.1, http://unesdoc.unesco.org/images/0018/001864/186462e.pdf (last access: 16 August 2021), 2009.

Fourteau, K., Hagenmuller, P., Roulle, J., and Domine, F.: On the use of heated needle probes for measuring snow thermal conductivity, Journal of Glaciology, 68, 705–719, https://doi.org/10.1017/jog.2021.127, 2022.

Fourteau, K., Freitag, J., Malinen, M., and Löwe, H.: Microstructure-based simulations of the viscous den-sification of snow and firn, EGUSphere [preprint], https://doi.org/10.5194/egusphere-2023-1928, 2023.
Gammon, P. H., Kiefte, H., Clouter, M. J., and Denner, W. W.: Elastic Constants of Artificial and Natural Ice Samples by Brillouin Spectroscopy, Journal of Glaciology, 29, 433–460, https://doi.org/10.3189/S0022143000030355, 1983.

Ganju, A., Satyawali, P. K., and Singh, A.: Snowcover Simulation Model: A Review, Defence Science Journal, 49, 419, https://unesdoc.unesco.org/ark:/48223/pf0000186462 (last access: 13 February 2026), 1999.

Gaume, J., Löwe, H., Tan, S., and Tsang, L.: Scaling laws for the mechanics of loose and cohesive granular materials based on Baxter's sticky hard spheres, Physical Review E, 96, 032914, https://doi.org/10.1103/PhysRevE.96.032914, 2017.

Gelebart, L., Derouillat, J., Doucet, N., Ouaki, F., Marano, A., and Duverge, J.: Amitex_FFTP, https://amitexfftp.github.io/AMITEX/general/index.html (last access: 13 February 2026), 2020.

Glen, J. W.: The Creep of Polycrystalline Ice, Proceedings of the Royal Society of London Series A, 228, 519–538, https://doi.org/10.1098/rspa.1955.0066, 1955.

Goldsby, D. L. and Kohlstedt, D. L.: Grain boundary sliding in fine-grained Ice I, Scripta Materialia, 37, 1399–1406, https://doi.org/10.1016/S1359-6462(97)00246-7, 1997.

Goldsby, D. L. and Kohlstedt, D. L.: Superplastic deformation of ice: Experimental observations, Journal of Geophysical Research: Solid Earth, 106, 11017–11030, https://doi.org/10.1029/2000JB900336, 2001.

Gorynina, O., Bartelt, P., and Gorynin, G.: One-Dimensional Visco-Elastic Finite Element Modeling of the Snow Creep, SSRN, https://doi.org/10.2139/ssrn.4781843, 2024.

Granger, R., Flin, F., Ludwig, W., Hammad, I., and Geindreau, C.: Orientation selective grain sublimation–deposition in snow under temperature gradient metamorphism observed with diffraction contrast tomography, The Cryosphere, 15, 4381–4398, https://doi.org/10.5194/tc-15-4381-2021, 2021.

Gray, J. M. N. T. and Morland, L. W.: The compaction of polar snow packs, Cold Regions Science and Technology, 23, 109–119, https://doi.org/10.1016/0165-232X(94)00010-U, 1995.

Green, D. J.: An introduction to the mechanical properties of ceramics, p. 348, Editor Cambridge University Press, ISBN 052159913X, 9780521599139, 336 pp., 1998.

Hagenmuller, P.: Modélisation du comportement mécanique de la neige à partir d'images microtomographiques, phdthesis, Université de Grenoble, https://theses.hal.science/tel-01230595 (last access: 13 February 226), 2014.

Hagenmuller, P. and Carmagnola, C.: Granulométrie de la neige de culture, Neige et Avalanches, 167, https://ens.syracuse.cloud/Default/doc/SYRACUSE/70156/granulometrie-de-la-neige-de-culture (last access: 13 February 226), 2020.

Hagenmuller, P., Calonne, N., Chambon, G., Flin, F., Geindreau, C., and Naaim, M.: Characterization of the snow microstructural bonding system through the minimum cut density, Cold Regions Science and Technology, 108, 72–79, https://doi.org/10.1016/j.coldregions.2014.09.002, 2014a.

Hagenmuller, P., Chambon, G., Flin, F., Morin, S., and Naaim, M.: Snow as a granular material: assessment of a new grain segmentation algorithm, Granular Matter, 16, 421–432, https://doi.org/10.1007/s10035-014-0503-7, 2014b.

Hagenmuller, P., Matzl, M., Chambon, G., and Schneebeli, M.: Sensitivity of snow density and specific surface area measured by microtomography to different image processing algorithms, The Cryosphere, 10, 1039–1054, https://doi.org/10.5194/tc-10-1039-2016, 2016.

Helfer, T., Michel, B., Proix, J.-M., Salvo, M., Sercombe, J., and Casella, M.: Introducing the open-source mfront code generator: Application to mechanical behaviours and material knowledge management within the PLEIADES fuel element modelling platform, Computers & Mathematics with Applications, 70, 994–1023, https://doi.org/10.1016/j.camwa.2015.06.027, 2015.

Hill, R.: Elastic properties of reinforced solids: Some theoretical principles, Journal of the Mechanics and Physics of Solids, 11, 357–372, https://doi.org/10.1016/0022-5096(63)90036-X, 1963.

Holman, L. E. and Leuenberger, H.: The relationship between solid fraction and mechanical properties of compacts – the percolation theory model approach, International Journal of Pharmaceutics, 46, 35–44, https://doi.org/10.1016/0378-5173(88)90007-5, 1988.

Huo, H., Chen, Q., Xiao, E., Li, H., Xu, H., Li, T., and Tang, X.: Long-Term One-Dimensional Compression Tests and Fractional Creep Model of Compacted Snow, Cold Regions Science and Technology, 228, 104326, https://doi.org/10.1016/j.coldregions.2024.104326, 2024.

Ignat, M. and Frost, H. J.: Grain boundary sliding in ice, Le Journal de Physique Colloques, 48, C1-189–C1-195, https://doi.org/10.1051/jphyscol:1987127, 1987.

Jacka, T. H.: The time and strain required for development of minimum strain rates in ice, Cold Regions Science and Technology, 8, 261–268, https://doi.org/10.1016/0165-232X(84)90057-0, 1984.

Johnson, J. B. and Hopkins, M. A.: Identifying microstructural deformation mechanisms in snow using discrete-element modeling, Journal of Glaciology, 51, 432–442, https://doi.org/10.3189/172756505781829188, 2005.

Keeler, C. M.: The Growth of Bonds and the Increase of Mechanical Strength in a Dry Seasonal Snow-Pack, Journal of Glaciology, 8, 441–450, https://doi.org/10.3189/S0022143000027027, 1969.

Kirchner, H. K., Michot, G., Narita, H., and Suzuki, T.: Snow as a foam of ice: Plasticity, fracture and the brittle-to-ductile transition, Philosophical Magazine A, 81, 2161–2181, https://doi.org/10.1080/01418610108217141, 2001.

Kojima, K.: Densification of seasonal snow cover, Physics of Snow and Ice, 1, 929–952, https://cir.nii.ac.jp/crid/1570854175146603264 (last access: 13 February 2026), 1967.

Kojima, K.: A field experiment on the rate of densification of natural snow layers under low stresses, in: Snow Mechanics Symposium, International Conference on Low Temperature Science, I. Conference on Physics of Snow and Ice, II. Conference on Cryobiology, 14–19 August 1966, Sapporo, Japan, https://cir.nii.ac.jp/crid/1570854175146603264 (last access: 13 February 2026), 1975.

Kominami, Y., Endo, Y., Niwano, S., and Ushioda, S.: Viscous compression model for estimating the depth of new snow, Annals of Glaciology, 26, 77–82, https://doi.org/10.3189/1998AoG26-1-77-82, 1998.

Lafaysse, M., Cluzet, B., Dumont, M., Lejeune, Y., Vionnet, V., and Morin, S.: A multiphysical ensemble system of nu-
merical snow modelling, The Cryosphere, 11, 1173–1198,
https://doi.org/10.5194/tc-11-1173-2017, 2017.

Lebensohn, R., Montagnat, M., Mansuy, P., Duval, P., Meyson-
nier, J., and Philip, A.: Modeling viscoplastic behavior
and heterogeneous intracrystalline deformation of colum-
nar ice polycrystals, Acta Materialia, 57, 1405–1415,
https://doi.org/10.1016/j.actamat.2008.10.057, 2009.

Lehning, M., Bartelt, P., Brown, B., Fierz, C., and Satyawali,
P.: A physical SNOWPACK model for the Swiss avalanche
warning: Part II, Snow microstructure, Cold Regions Science
and Technology, 35, 147–167, https://doi.org/10.1016/S0165-
232X(02)00073-3, 2002.

Lejeune, Y., Dumont, M., Panel, J.-M., Lafaysse, M., Lapalus, P., Le Gac, E., Lesaffre, B., and Morin, S.: 57 years (1960–2017) of snow and meteorological observations from a mid-altitude mountain site (Col de Porte, France, 1325 m of altitude), Earth System Science Data, 11, 71–88, https://doi.org/10.5194/essd-11-71-2019, 2019.

Louchet, F.: Dislocations and plasticity in
ice, Comptes Rendus Physique, 5, 687–698,
https://doi.org/10.1016/j.crhy.2004.09.001, 2004.

Lundin, J. M. D., Stevens, C. M., Arthern, R., Buizert, C., Orsi, A., Ligtenberg, S. R. M., Simonsen, S. B., Cummings, E., Essery, R., Leahy, W., Harris, P., Helsen, M. M., and Waddington, E. D.: Firn Model Intercomparison Experiment (FirnMICE), Journal of Glaciology, 63, 401–422, https://doi.org/10.1017/jog.2016.114, Press, 2017.

Lundy, C. C., Brown, R. L., Adams, E. E., Birkeland, K. W., and Lehning, M.: A statistical validation of the snowpack model in a Montana climate, Cold Regions Science and Technology, 33, 237–246, https://doi.org/10.1016/S0165-232X(01)00038-6, 2001.

Löwe, H., Zaiser, M., Mösinger, S., and Schleef, S.: Snow Me-
chanics Near the Ductile-Brittle Transition: Compressive Stick-
Slip and Snow Microquakes, Geophysical Research Letters,
47, e2019GL085491, https://doi.org/10.1029/2019GL085491, 2020.

McClung, D. M. and Larsen, J. O.: Snow creep pressures: Effects of structure boundary conditions and snowpack properties compared with field data, Cold Regions Science and Technology, 17, 33–47, https://doi.org/10.1016/S0165-232X(89)80014-X, 1989.

Mellor, M.: A review of basic snow mechanics, IAHS Publ., 114, 251–291, 1974.

Mellor, M. and Cole, D. M.: Deformation and failure of ice under constant stress or constant strain-rate, Cold Regions Science and Technology, 5, 201–219, https://doi.org/10.1016/0165-232X(82)90015-5, 1982.

Mellor, M. and Testa, R.: Creep of Ice under Low Stress, Journal of Glaciology, 8, 147–152, https://doi.org/10.3189/S0022143000020815, 1969a.

Mellor, M. and Testa, R.: Effect of Temperature on the Creep of Ice, Journal of Glaciology, 8, 131–145, https://doi.org/10.3189/S0022143000020803, 1969b.

Meyssonnier, J., Philip, A., Capolo, L., and Mansuy, P.: Experimen-
tal studies of the viscoplasticity of ice and snow, in: Mechanics
of Natural Solids, edited by: Kolymbas, D. and Viggiani, G.,
203–221, Springer, Berlin, Heidelberg, ISBN 978-3-642-03578-
4, https://doi.org/10.1007/978-3-642-03578-4_9, 2009.

Montagnat, M., Castelnau, O., Bons, P. D., Faria, S. H., Gagliar-
dini, O., Gillet-Chaulet, F., Grennerat, F., Griera, A., Leben-

sohn, R. A., and Moulinec, H.: Multiscale modeling of ice de-
formation behavior, Journal of Structural Geology, 61, 78–108,
https://doi.org/10.1016/j.jsg.2013.05.002, 2014.

Moos, M. v., Bartelt, P., Zweidler, A., and Bleiker, E.: Triaxial tests on snow at low strain rate. Part I. Experimental device, Journal of Glaciology, 49, 81–90, https://doi.org/10.3189/172756503781830881, 2003.

Morin, S., Horton, S., Techel, F., Bavay, M., Coléou, C., Fierz, C., Gobiet, A., Hagenmuller, P., Lafaysse, M., Ližar, M., Mitrerer, C., Monti, F., Müller, K., Olefs, M., Snook, J. S., van Herwijnen, A., and Vionnet, V.: Application of physical snowpack models in support of operational avalanche hazard forecasting: A status report on current implementations and prospects for the future, Cold Regions Science and Technology, 170, 102910, https://doi.org/10.1016/j.coldregions.2019.102910, 2020.

Morris, E. M. and Wingham, D. J.: Densification of polar snow: Measurements, modeling, and implications for altimetry, Journal of Geophysical Research: Earth Surface, 119, 349–365, https://doi.org/10.1002/2013JF002898, 2014.

Narita, H.: Mechanical behaviour and structure of snow under uniaxial tensile stress, Journal of Glaciology, 26, 275–282, https://doi.org/10.3189/S0022143000010819, 1980.

Narita, H.: An experimental study on tensile fracture of snow, Con-
tributions from the institute of Low Temperature Science, Insti-
tute of Low Temperature Science, Hokkaido University, 32, 1–
37, https://eprints.lib.hokudai.ac.jp/repo/huscap/all/20246/A32_
p1-37.pdf (last access: 13 February 2026), 1984.

Navarre, J. P.: Modele unidimensionnel d'evolution de la neige deposee, Société météorologique de France, La Météorologie, Viéme série, 12/75, 1975.

O'Connor, D. and Haehnel, R.: A generalized approach for modeling creep of snow foundations, Tech. rep., Engineer Research and Development Center (U.S.), https://doi.org/10.21079/11681/36593, 2020.

Peinke, I., Hagenmuller, P., Andò, E., Chambon, G., Flin, F., and Roulle, J.: Experimental Study of Cone Penetration in Snow Using X-Ray Tomography, Frontiers in Earth Science, 8, https://doi.org/10.3389/feart.2020.00063, 2020.

Proksch, M., Löwe, H., and Schneebeli, M.: Density, specific sur-
face area, and correlation length of snow measured by high-
resolution penetrometry, Journal of Geophysical Research: Earth
Surface, 120, 346–362, https://doi.org/10.1002/2014JF003266,
2015.

Quéno, L., Vionnet, V., Dombrowski-Etchevers, I., Lafaysse, M., Dumont, M., and Karbou, F.: Snowpack modelling in the Pyrenees driven by kilometric-resolution meteorological forecasts, The Cryosphere, 10, 1571–1589, https://doi.org/10.5194/tc-10-1571-2016, 2016.

Ramseier, R. O.: Growth and mechanical properties of river and lake ice, PhD Thesis, Laval University, Canada, 1972.

Ritter, J., Löwe, H., and Gaume, J.: Microstructural controls of anticrack nucleation in highly porous brittle solids, Scientific Reports, 10, 12383, https://doi.org/10.1038/s41598-020-67926-2, 2020.

Roscoat, S. R. D., King, A., Philip, A., Reischig, P., Ludwig, W., Flin, F., and Meyssonnier, J.: Analysis of Snow Microstructure by Means of X-Ray Diffraction Contrast Tomography, Advanced Engineering Materials, 13, 128–135, https://doi.org/10.1002/adem.201000221, 2011.
Royer, A., Picard, G., Vargel, C., Langlois, A., Gouttevin, I., and Dumont, M.: Improved Simulation of Arctic Circumpolar Land Area Snow Properties and Soil Temperatures, Frontiers in Earth Science, 9, https://doi.org/10.3389/feart.2021.685140, 2021.

Scapozza, C. and Bartelt, P.: Triaxial tests on snow at low strain rate, Part II. Constitutive behaviour, Journal of Glaciology, 49, 91–101, https://doi.org/10.3189/172756503781830890, 2003a.

Scapozza, C. and Bartelt, P. A.: The influence of temperature on the small-strain viscous deformation mechanics of snow: a comparison with polycrystalline ice, Annals of Glaciology, 37, 90–96, https://doi.org/10.3189/172756403781815410, 2003b.

Schleef, S., Löwe, H., and Schneebeli, M.: Hot-pressure sintering of low-density snow analyzed by X-ray microtomography and in situ microcompression, Acta Materialia, 71, 185–194, https://doi.org/10.1016/j.actamat.2014.03.004, 2014a.

Schleef, S., Löwe, H., and Schneebeli, M.: Influence of stress, temperature and crystal morphology on isothermal densification and specific surface area decrease of new snow, The Cryosphere, 8, 1825–1838, https://doi.org/10.5194/tc-8-1825-2014, 2014b.

Schneider, M., Merkert, D., and Kabel, M.: FFT-based homogenization for microstructures discretized by linear hexahedral elements, International Journal for Numerical Methods in Engineering, 109, 1461–1489, https://doi.org/10.1002/nme.5336, 2017.

Schulson, E. M. and Duval, P.: Creep and Fracture of Ice, Cambridge University Press, ISBN 978-0-521-80620-6, 2009.

Schöttner, J., Walet, M., Rosendahl, P., Weissgraeber, P., Adam, V., Walter, B., Rheinschmidt, F., Löwe, H., Schweizer, J., and Herwijnen, A. V.: On the compressive strength of weak snow layers of depth hoar, Journal of Glaciology, 71, e54, https://doi.org/10.1017/jog.2025.16, 2025.

Shinojima, K.: Study on the visco-elastic deformation of deposited snow, Physics of snow and ice, The Institute of Low Temperature Science, Hokkaido University, https://cir.nii.ac.jp/crid/1572824500219819520 (last access: 13 February 2026), 1967.

Simson, A., Löwe, H., and Kowalski, J.: Elements of future snow-pack modeling – Part 2: A modular and extendable Eulerian–Lagrangian numerical scheme for coupled transport, phase changes and settling processes, The Cryosphere, 15, 5423–5445, https://doi.org/10.5194/tc-15-5423-2021, 2021.

Srivastava, P. K., Chandel, C., Mahajan, P., and Pankaj, P.: Prediction of anisotropic elastic properties of snow from its microstructure, Cold Regions Science and Technology, 125, 85–100, https://doi.org/10.1016/j.coldregions.2016.02.002, 2016.

Stauffer, D. and Aharony, A.: Introduction To Percolation Theory: Second Edition, Taylor & Francis, London, 2nd Edn., ISBN 978-1-315-27438-6, https://doi.org/10.1201/9781315274386, 2018.

Steinkogler, W., Fierz, C., Lehning, M., and Obleitner, F.: Systematic Assessment of New Snow Settlement in Snowpack, International Snow Science Workshop, Davos 2009, Proceedings, 132–135, https://arc.lib.montana.edu/snow-science/item/211 (last access: 13 February 2026), 2009.

Stoffel, M. and Bartelt, P.: Modelling Snow Slab Release Using a Temperature-Dependent Viscoelastic Finite Element Model with Weak Layers, Surveys in Geophysics, 24, 417–430, https://doi.org/10.1023/B:GEOP.0000006074.56474.43, 2003.

Sturm, M. and Holmgren, J.: Differences in compaction behavior of three climate classes of snow, Annals of Glaciology, 26, 125–130, https://doi.org/10.3189/1998AoG26-1-125-130, 1998.

Sundu, K., Freitag, J., Fourteau, K., and Löwe, H.: A microstructure-based parameterization of the effective anisotropic elasticity tensor of snow, firn, and bubbly ice, The Cryosphere, 18, 1579–1596, https://doi.org/10.5194/tc-18-1579-2024, 2024a.

Sundu, K., Ottersberg, R., Jaggi, M., and Löwe, H.: A grain-size driven transition in the deformation mechanism in slow snow compression, Acta Materialia, 262, 119359, https://doi.org/10.1016/j.actamat.2023.119359, 2024b.

Suquet, P., Moulinec, H., Castelnau, O., Montagnat, M., Lallelec, N., Grennerat, F., Duval, P., and Brenner, R.: Multiscale modeling of the mechanical behavior of polycrystalline ice under transient creep, Procedia IUTAM, 3, 76–90, https://doi.org/10.1016/j.piutam.2012.03.006, 2012.

Teufelsbauer, H.: A two-dimensional snow creep model for alpine terrain, Natural Hazards, 56, 481–497, https://doi.org/10.1007/s11069-010-9515-8, 2011.

Theile, T., Löwe, H., Theile, T. C., and Schneebeli, M.: Simulating creep of snow based on microstructure and the anisotropic deformation of ice, Acta Materialia, 59, 7104–7113, https://doi.org/10.1016/j.actamat.2011.07.065, 2011.

Touzeau, A., Landais, A., Morin, S., Arnaud, L., and Picard, G.: Numerical experiments on vapor diffusion in polar snow and firn and its impact on isotopes using the multi-layer energy balance model Crocus in SURFEX v8.0, Geoscientific Model Development, 11, 2393–2418, https://doi.org/10.5194/gmd-11-2393-2018, 2018.

Treverrow, A., Budd, W. F., Jacka, T. H., and Warner, R. C.: The tertiary creep of polycrystalline ice: experimental evidence for stress-dependent levels of strain-rate enhancement, Journal of Glaciology, 58, 301–314, https://doi.org/10.3189/2012JoG11J149, 2012.

Underwood, E. E.: Quantitative Stereology, Addison-Wesley series in metallurgy and materials, ISSN 0515-3972, Addison-Wesley Publishing Company, l'Université de Californie, ISBN 0201076500, 9780201076509, 274 pp., 1970.

van Kampenhout, L., Lenaerts, J. T. M., Lipscomb, W. H., Sacks, W. J., Lawrence, D. M., Slater, A. G., and van den Broeke, M. R.: Improving the Representation of Polar Snow and Firn in the Community Earth System Model, Journal of Advances in Modeling Earth Systems, 9, 2583–2600, https://doi.org/10.1002/2017MS000988, 2017.

Verjans, V., Leeson, A. A., Stevens, C. M., MacFerrin, M., Noël, B., and van den Broeke, M. R.: Development of physically based liquid water schemes for Greenland firn-densification models, The Cryosphere, 13, 1819–1842, https://doi.org/10.5194/tc-13-1819-2019, 2019.

Vionnet, V., Brun, E., Morin, S., Boone, A., Faroux, S., Le Moigne, P., Martin, E., and Willemet, J.-M.: The detailed snowpack scheme Crocus and its implementation in SURFEX v7.2, Geoscientific Model Development, 5, 773–791, https://doi.org/10.5194/gmd-5-773-2012, 2012.

Védrine, L. and Hagenmuller, P.: Supplementary data for “Revisiting snow settlement with microstructural knowledge”, Zenodo [data set], https://doi.org/10.5281/zenodo.16778695, 2025.

Védrine, L., Hagenmuller, P., Gélébart, L., Montagnat, M., and Bernard, A.: Role of Ice Mechanics on Snow Viscoplasticity, Geophysical Research Letters, 51, e2023GL107676, https://doi.org/10.1029/2023GL107676, 2024.
Védrine, L., Hagenmuller, P., Gélébart, L., Montagnat, M., and Löwe, H.: Sensitivity of the viscoplasticity of polycrystals to porosity and pore-to-crystal size ratio, Acta Materialia, 301, 121507, https://doi.org/10.1016/j.actamat.2025.121507, 2025.

Wautier, A., Geindreau, C., and Flin, F.: 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, 8031–8041, https://doi.org/10.1002/2015GL065227, 2015.

Wautier, A., Geindreau, C., and Flin, F.: Numerical homogenization of the viscoplastic behavior of snow based on X-ray tomography images, The Cryosphere, 11, 1465–1485, https://doi.org/10.5194/tc-11-1465-2017, 2017.

Weikusat, I., Kuiper, E.-J. N., Pennock, G. M., Kipfstuhl, S., and Drury, M. R.: EBSD analysis of subgrain boundaries and dislocation slip systems in Antarctic and Greenland ice, Solid Earth, 8, 883–898, https://doi.org/10.5194/se-8-883-2017, 2017.

Wever, N., Schmid, L., Heilig, A., Eisen, O., Fierz, C., and Lehning, M.: Verification of the multi-layer SNOWPACK model with different water transport schemes, The Cryosphere, 9, 2271–2293, https://doi.org/10.5194/tc-9-2271-2015, 2015.

Willot, F.: Fourier-based schemes for computing the mechanical response of composites with accurate local fields, Comptes Rendus Mécanique, 343, 232–245, https://doi.org/10.1016/j.crme.2014.12.005, 2015.

Woolley, G. J., Rutter, N., Wake, L., Vionnet, V., Derksen, C., Essery, R., Marsh, P., Tutton, R., Walker, B., Lafaysse, M., and Pritchard, D.: Multi-physics ensemble modelling of Arctic tundra snowpack properties, The Cryosphere, 18, 5685–5711, https://doi.org/10.5194/tc-18-5685-2024, 2024.

Yosida, Z.: Physical studies on deposited snow thermal properties, Contributions from the Institute of Low Temperature Science, Ser. A., 7, 19–74, https://cir.nii.ac.jp/crid/1570009749325699456 (last access: 13 February 2026), 1955.