Revisiting snow settlement with microstructural knowledge

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 snowpack evolution (Brun et al., 1992; Lehning et al., 2002; Simson 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 paleoclimatology (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 $\dot{\mathit{ϵ}}$ (Glen, 1955):

$$\begin{array}{}\text{(1)}& \dot{\mathit{ϵ}}={\dot{\mathit{ϵ}}}_{\mathrm{0}}{\left({\displaystyle \frac{\mathit{\sigma }}{{\mathit{\sigma }}_{\mathrm{0}}}}\right)}^n\end{array}$$

with ${\dot{\mathit{ϵ}}}_{\mathrm{0}}=\mathrm{1}$ s⁻¹, n the stress exponent and σ₀ 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 η linearly relates the two directly-measurable macroscopic quantities, i.e., stress σ and strain rate $\dot{\mathit{ϵ}}$:

$$\begin{array}{}\text{(2)}& \dot{\mathit{ϵ}}={\displaystyle \frac{\mathit{\sigma }}{\mathit{\eta }}}.\end{array}$$

This model implicitly assumes that n=1 in Eq. (1) and corresponds to a Newtonian compressible viscous fluid. Several authors have encapsulated into η 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 η does not depend on σ and $\dot{\mathit{ϵ}}$, 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=\mathrm{4}\pm \mathrm{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≈ 1.8 for low-density snow, and regime II with n≈ 3.5 for higher density snow, with a transition density depending on temperature. In line with this finding, Schleef et al. (2014a) reported n≈ 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 $-{\mathrm{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≈1.9 for small grains to n≈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 σ₀) 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 intracrystalline 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∼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 η (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 intra-crystalline 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 σ₀ (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 Carmagnola, 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, $\mathrm{3}/\left(\mathrm{SSA}{\mathit{\rho }}_{\mathrm{ice}}\right)$, remains between 6 and 12 voxels. A cubic sub-volume of $\mathrm{250}\times \mathrm{250}\times \mathrm{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.

Hagenmuller et al. (2016)Peinke et al. (2020)Hagenmuller and Carmagnola (2020)Fourteau et al. (2022)Bernard et al. (2023)

Table 1 Sample origin and reference papers, type and number of samples and voxels size (range) of images used in this study.

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 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 (RGlr), Faceted Crystals (FC), Depth Hoar (DH), Machine Made (MM). The microstructures enclosed by a black circle correspond to those shown in Fig. 2.

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.

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{\mathit{\gamma }}}_{\mathrm{0}}^{\left(k\right)}$, a stress exponent n^(k), and a critical shear stress ${\mathit{\tau }}_{\mathrm{0}}^{\left(k\right)}$. Different parametrizations of the slip systems have been proposed in the literature for different 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 (${\mathit{\sigma }}_{\mathrm{0}}(\mathrm{\Phi }=\mathrm{1})=\mathrm{272}\mathrm{MPa}$ and n_suquet=2.94) as proposed by Budd and Jacka (1989) and Castelnau et al. (1996) at T₀=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{\mathit{\gamma }}}_{\mathrm{0}}^{\left(k\right)}\left(T\right)={\dot{\mathit{\gamma }}}_{\mathrm{0}}^{\left(k\right)}\left(T_{\mathrm{0}}\right)A\left(T\right)$ with

$$\begin{array}{}\text{(3)}& A\left(T\right)=\exp \left[-{\displaystyle \frac{Q}{R}}\left({\displaystyle \frac{\mathrm{1}}{T}}-{\displaystyle \frac{\mathrm{1}}{T_{\mathrm{0}}}}\right)\right],\end{array}$$

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

Figure 3 Numerical boundary conditions: uniaxial compression test with imposed strain rate.

For the boundary conditions, we prescribed vertical uniaxial compression with periodic boundary conditions (Fig. 3). We imposed a positive average vertical strain rate ${\dot{\mathit{\epsilon }}}_{zz}>\mathrm{0}$ and zero average lateral stresses, ${\mathit{\sigma }}_{xx}={\mathit{\sigma }}_{yy}=\mathrm{0}$, and zero average shear stresses, ${\mathit{\sigma }}_{xy}={\mathit{\sigma }}_{yz}={\mathit{\sigma }}_{xz}=\mathrm{0}$. Here, the macroscopic strain and stress tensors are defined as volume average on the sample: $\mathit{\epsilon }=\frac{\mathrm{1}}{\left|V\right|}{\int }_V\mathit{\epsilon }\left(x\right)\mathrm{d}V$ and $\mathit{\sigma }=\frac{\mathrm{1}}{\left|V\right|}{\int }_V\mathit{\sigma }\left(x\right)\mathrm{d}V$, 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 σ_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 $\mathrm{250}\times \mathrm{250}\times \mathrm{250}$ voxels, on a high-performance computer with 256 IntelXeon Gold 6130 cores). However, the solver benefits from the MPI implementation.

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⁻⁷ to ${\mathrm{10}}^{-\mathrm{5}}{\text{s}}^{-\mathrm{1}}$. We observe three different regimes of the stress-strain response (Fig. 4). First, the axial stress σ_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, σ_zz deviates from the elasticity, and the viscoplasticity is progressively activated. After this transient regime, σ_zz eventually reaches a constant value, the flow stress σ_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 typical 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 σ_Y.

Figure 4 Evolution of the vertical stress σ_zz as the function of the strain ϵ_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 Φ=0.48, and is shown in Fig. 3.

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 σ₀. 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:

$$\begin{array}{}\text{(4)}& \dot{\mathit{\epsilon }}={\dot{\mathit{\epsilon }}}_{\mathrm{0}}\left(T_{\mathrm{0}}\right)A\left(T\right){\left({\displaystyle \frac{\mathit{\sigma }}{{\mathit{\sigma }}_{\mathrm{0}}}}\right)}^n\end{array}$$

with n and σ₀ 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 $\mathrm{4}\times {\mathrm{10}}^{-\mathrm{4}}$ s⁻¹ 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{\mathit{ϵ}},\mathit{\sigma })$ and essential information such as initial ice volume fraction Φ, 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).

Scapozza and Bartelt (2003b)Theile et al. (2011)Schleef et al. (2014b)Bernard (2023)Bernard (2023)Bernard et al. (2023)

Table 2 Characteristic of the experimental data of snow creep used.

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.

It is impossible to derive two parameters, namely the reference stress σ₀ and stress exponent n, from an individual measurement point $(\dot{\mathit{ϵ}},\mathit{\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 σ₀(n) as:

$$\begin{array}{}\text{(5)}& {\mathit{\sigma }}_{\mathrm{0}}\left(n\right)=\mathit{\sigma }\times {\left({\displaystyle \frac{\dot{\mathit{ϵ}}}{{\dot{\mathit{\epsilon }}}_{\mathrm{0}}\left(T_{\mathrm{0}}\right)A\left(T\right)}}\right)}^{-\frac{\mathrm{1}}{n}}.\end{array}$$

If the power-law is well-formed, σ₀ 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 σ₀(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 Φ 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 ${\mathit{\sigma }}_{\mathrm{0}}={\mathit{\sigma }}_{\mathrm{0}}(n,\mathrm{\Phi })$ and n=n(Φ). Consequently, the exponent n(Φ) is the one that minimizes the scatter σ₀(n,Φ) for a given narrow range of Φ 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 ℒ(model parameters) defined as:

$$\begin{array}{}\text{(6)}& \mathcal{L}\left(\text{model parameters}\right)={\displaystyle \frac{\mathrm{1}}{N_{\text{dataset}}}}\sum _{i\in \text{dataset}}{\displaystyle \frac{{\left(\stackrel{\mathrm{̃}}{{\mathit{\sigma }}_{\mathrm{0}}}\left({\mathrm{\Phi }}_i\right)-{\mathit{\sigma }}_{{\mathrm{0}}_i}\right)}^{\mathrm{2}}}{{\mathit{\sigma }}_{{\mathrm{0}}_i}^{\mathrm{2}}}}\end{array}$$

and where ${\mathit{\sigma }}_{{\mathrm{0}}_i}$ and $\stackrel{\mathrm{̃}}{{\mathit{\sigma }}_{\mathrm{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 σ is expressed as:

$$\begin{array}{}\text{(7)}& {\displaystyle \frac{\mathrm{d}D}{D}}=-{\displaystyle \frac{\mathit{\sigma }}{\mathit{\eta }}}\mathrm{d}t\end{array}$$

where dt is the model time step and η (in MPa s) denotes the compaction viscosity. The compaction viscosity is typically parametrised as a function of temperature T (in °C) and density ρ (in kg m⁻³).

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):

$$\begin{array}{}\text{(8)}& {\mathit{\eta }}_{\mathrm{BR}\mathrm{92}}=\mathrm{4}{\mathit{\eta }}_{\mathrm{0}}{\displaystyle \frac{\mathit{\rho }}{c_{\mathit{\eta }}}}\exp (-a_{\mathit{\eta }}T+b_{\mathit{\eta }}\mathit{\rho })\end{array}$$

where η₀=7.62237 MPa, a_η=0.1 °C⁻¹, b_η=0.023 m³ kg⁻¹, and c_η=250 kg m⁻³.

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):

$$\begin{array}{}\text{(9)}& {\mathit{\eta }}_{\mathrm{T}\mathrm{11}}=\mathrm{5}\times {\mathrm{10}}^{-\mathrm{8}}{\mathit{\rho }}^{(-\mathrm{0.0371}T+\mathrm{4.4})}\left({\mathrm{10}}^{-\mathrm{4}}\exp \left(\mathrm{0.018}\mathit{\rho }\right)+\mathrm{1}\right)\end{array}$$

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:

$$\begin{array}{}\text{(10)}& {\mathit{\eta }}_{\mathrm{K}\mathrm{75}}=\mathrm{7}\times {\mathrm{10}}^{-\mathrm{9}}{\mathit{\rho }}^{(\mathrm{4.75}-T/\mathrm{40})}\end{array}$$

Although viscosity is a material parameter that characterises deformation in a specific direction, it is commonly assumed that for densities below 450 kg m⁻³, 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.

Figure 5 Parameters of the homogenized model: (a) stress exponent n and (b) reference stress σ₀ as a function of the solid fraction Φ.

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 Φ. The stress exponent ranges between 2.1 and 2.2 for the considered microstructures with Φ between 0.1 and 0.6. These values are lower than the exponent of dense ice ($n_{{\mathrm{ice}}_{\mathrm{suquet}}}=\mathrm{2.94}$) and are close to the basal slip exponent (n_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 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 inter-crystalline 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:

$$\begin{array}{}\text{(11)}& {\mathit{\sigma }}_{\mathrm{0}}={\mathit{\sigma }}_{\mathrm{0}}(\mathrm{\Phi }=\mathrm{1}){\mathrm{\Phi }}_{\mathrm{r}}^m.\end{array}$$

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

$$\begin{array}{}\text{(12)}& {\mathrm{\Phi }}_{\mathrm{r}}={\displaystyle \frac{\mathrm{\Phi }-{\mathrm{\Phi }}_t}{\mathrm{1}-{\mathrm{\Phi }}_t}}\end{array}$$

The quantity ${\mathit{\sigma }}_{\mathrm{0}}(\mathrm{\Phi }=\mathrm{1})=\mathrm{272}\mathrm{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 ℒ=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 ϕ^−nm (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 (ℒ=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 Φ: σ₀∼Φ^3.96. 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 ℒ(n) (defined as the relative standard deviation, see Eq. 6) of σ₀ around a power-law of Φ of the form $\stackrel{\mathrm{̃}}{{\mathit{\sigma }}_{\mathrm{0}}}\left(n\right)={\mathit{\sigma }}_{\mathrm{0}}(\mathrm{\Phi }=\mathrm{1}){\mathrm{\Phi }}_{\mathrm{r}}^m$, with Φ_r the rescaled solid fraction (Eq. 12) and, m and Φ_t two optimal fitting parameters depending on n.

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.

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 σ₀ (ℒ=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 ℒ=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_µCT=2.15, the dispersion is reduced to ℒ=0.14, with optimized parameters Φ_t=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 ℒ as a function of n. ℒ decreases when n increases from 1 to 2.35. For n in [2.35, 2.59], ℒ remains nearly constant with a weak minimum at n=2.53. For n>2.59, ℒ 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, $\mathcal{L}(n=\mathrm{2.15})=\mathrm{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 (ℒ=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 Φ^m 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 ([$\sim -\mathrm{18.5},\sim -\mathrm{11},\sim -\mathrm{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 deformation mechanisms may take place (Mellor and Testa, 1969b).

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.

Figure 8 Normalized residual $R=\frac{{\dot{\mathit{ϵ}}}_{V\mathrm{25}}-{\dot{\mathit{ϵ}}}_{\mathrm{Br}\mathrm{92}}}{{\dot{\mathit{ϵ}}}_{\mathrm{Br}\mathrm{92}}+{\dot{\mathit{ϵ}}}_{\mathrm{V}\mathrm{25}}}$ between the V25 and Br92 model for varying solid fractions and stresses. The stress is converted into its snow water equivalent.

Figure 7 illustrates the reference stress σ₀ 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 around 0.4, suggesting a potential evolution of the morphology parameter m, which is here assumed constant for Φ in the range [Φ_t, 1]. Indeed, Coble (1961) showed that for Φ 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):

$$\begin{array}{}\text{(13)}& \dot{\mathit{ϵ}}(\mathit{\sigma },T)=e^{\frac{-Q}{R}\left(\frac{\mathrm{1}}{T}-\frac{\mathrm{1}}{T_{\mathrm{0}}}\right)}{\left({\displaystyle \frac{\mathit{\sigma }}{{\mathit{\sigma }}_{\mathrm{0}(\mathrm{\Phi }=\mathrm{1})}{\mathrm{\Phi }}_{\mathrm{r}}^m}}\right)}^n,\end{array}$$

with $Q=\mathrm{69.1}\mathrm{kJ}{\mathrm{mol}}^{-\mathrm{1}}$, T₀=263 K, n=2.15, m=3.195 and Φ_r the rescaled solid fraction (Eq. 12) with Φ_t=0.025. The evaluation is focused on stress - solid fraction (σ,Φ) 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 σ is related to strain rate as $\dot{\mathit{ϵ}}=\dot{h}/h$. Note that the overburden stress is defined as the gravitational acceleration (g=9.81 m s⁻²) 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 η on solid fraction Φ. However, substantial discrepancies arise when moving beyond these “standard” conditions: Br92 predicts higher strain rates than V25 for Φ>0.08 and stresses smaller than those observed at Col de Porte. Conversely, for very low density samples (Φ<0.08) and for stresses higher than the “standard” ones, Br92 predicts a lower strain rate compared to V25.

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.

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.

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.

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 kg m⁻³ (Φ=0.076, σ=0.17 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 underestimated 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 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 intercrystalline 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, Φ≤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{\mathit{ϵ}}={\dot{\mathit{ϵ}}}_{\mathrm{0}}(\mathit{\sigma }/{\mathit{\sigma }}_{\mathrm{0}}{)}^n$. The stress exponent n ranged from 2.1 to 2.2. The reference stress σ₀ was found to be a function of the solid fraction Φ, following the relationship ${\mathit{\sigma }}_{\mathrm{0}}\sim {\mathrm{\Phi }}_{\mathrm{r}}^m$. 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 ($\mathit{\sigma },\dot{\mathit{ϵ}}$) 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 η, in snowpack models, implicitly accounts for both the material viscosity variation with Φ and the site-specific correlation between the stress and Φ. As a result, η 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 ε is decomposed into an elastic part (ε^e) and a viscoplastic part (ε^vp), such that: $\mathit{\epsilon }={\mathit{\epsilon }}^e+{\mathit{\epsilon }}^{vp}$.

The elastic strain ε^e is related to the stress tensor σ and the fourth-order stiffness tensor C by the relation:

$$\begin{array}{}\text{(A1)}& \mathit{\sigma }=C:{\mathit{\epsilon }}^e.\end{array}$$

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₁₁=13.9 GPa, C₃₃=15.0 GPa, C₄₄=3.0 GPa, C₁₂=7.1 GPa, C₁₃=5.8 GPa.

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

$$\begin{array}{}\text{(A2)}& {\dot{\mathit{\epsilon }}}^{vp}=\sum _{k=\mathrm{1}}^N{\dot{\mathit{\gamma }}}_{\mathrm{0}}^{\left(k\right)}\left(T\right){\left({\displaystyle \frac{\left|{\mathit{\tau }}^{\left(k\right)}\right|}{{\mathit{\tau }}_{\mathrm{0}}^{\left(k\right)}}}\right)}^{n^{\left(k\right)}}\mathrm{sgn}\left({\mathit{\tau }}^{\left(k\right)}\right){\mathit{\mu }}^{\left(k\right)}\end{array}$$

where n^(k) is the stress exponent, ${\dot{\mathit{\gamma }}}_{\mathrm{0}}^{\left(k\right)}$ is the reference shear rate, ${\mathit{\tau }}_{\mathrm{0}}^{\left(k\right)}$ is the critical shear stress for slip system (k), and ${\mathit{\mu }}^{\left(k\right)}=n^{\left(k\right)}{\otimes }_{\mathrm{s}}{b}^{\left(k\right)}$ is the Schmid tensor defined by the (symmetric) dyadic product (⊗_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 ${\mathit{\tau }}_{\mathrm{0}}^{\left(k\right)}$ are assumed to be constant with time (no hardening), and the reference shear rate is set to ${\dot{\mathit{\gamma }}}_{\mathrm{0}}^{\left(k\right)}=\mathrm{1}$ s⁻¹. 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=\mathrm{7.8}\times {\mathrm{10}}^{-\mathrm{8}}$ MPa⁻¹ s⁻³ and n_poly=2.94 the parameters of the Glen law for polycrystalline ice $\dot{\mathit{ϵ}}=A{\mathit{\sigma }}^{n_{\mathrm{poly}}}$). The parametrisations of the crystal plasticity model are summarised in Table A1.

Table A1 Single crystal ice parameters of crystal plasticity models at −10 °C (Suquet et al., 2012).

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 ϕ=0.4, the snow microstructures exhibit a surface ratio of approximately β≈0.1, whereas the Voronoi microstructures exhibit a value of about β≈0.35.

Figure B1 Contiguity β and surface ratio r as the function of solid fraction Φ. The contiguity Underwood (1970) is defined as $\mathit{\beta }={\displaystyle \frac{r}{\mathrm{1}+r}}$.

Appendix C: Evolution of ℒ as a function of stress exponent n.

Figure C1 Evolution of the dispersion ℒ as a function of stress exponent n at a temperature of −10 °C.