Role of Ice Mechanics on Snow Viscoplasticity

The porous structure of snow becomes denser with time under gravity, primarily due to the creep of its ice matrix with viscoplasticity. Despite investigation of this behavior at the macroscopic scale, the driving microscopic mechanisms are still not well understood. Thanks to high‐performance computing and dedicated solvers, we modeled snow elasto‐viscoplasticity with 3D images of its microstructure and different mechanical models of ice. The comparison of our numerical experiments to oedometric compression tests measured by tomography showed that ice in snow rather behaves as a heterogeneous set of ice crystals than as homogeneous polycrystalline ice. Similarly to dense ice, the basal slip system contributed at most, in the simulations, to the total snow deformation. However, in the model, the deformation accommodation between crystals was permitted by the pore space and did not require any prismatic and pyramidal slips, whereas the latter are pre‐requisite for the simulation of dense ice.


Introduction
Once on the ground, snow naturally settles under its weight.Its density typically increases from 100 kg m 3 just after its deposition on the ground to 550 kg m 3 at the end of the winter season at mid-latitude regions, up to the density of ice, 917 kg m 3 , for buried perennial snow on glaciers or ice sheets (Alley, 1987;Arnaud, Lipenkov, et al., 1998).The primary driving mechanisms of this settlement are the deformation of the ice skeleton of snow and the subsequent pore volume reduction (e.g., Yosida et al., 1958).Understanding the viscoplastic behavior of snow is crucial to predict its densification, which is required to model the snowpack evolution (Lehning et al., 2002;Simson et al., 2021;Vionnet et al., 2012) or the pore close-off in ice cores (Gregory et al., 2014), which in turn relates to many applications such as avalanche forecasting (Morin et al., 2020), hydrology (DeBeer & Pomeroy, 2017) or paleo-climatology (Barnola et al., 1987).
The snow microstructure comprises open pores and sintered crystals of hexagonal ice (Ih).From a mechanical standpoint, the contribution of air is negligible, and snow mechanics thus depends solely on the spatial arrangement of these crystals and the deformation mechanisms they undergo.Under an overburden, isolated ice crystals essentially deform plastically by gliding of basal dislocations (Duval et al., 1983b).Non-basal dislocations were introduced in models to accommodate basal glide, but are poorly characterized (Hondoh, 2000).Their creep rate is at least two orders of magnitude lower than the one of the basal glide (Duval et al., 1983b).Ice crystals thus exhibit an extremely anisotropic viscoplastic behavior.When the crystals are sintered together, the grain boundaries act as obstacles for dislocation movement, which results in hardening and more complex deformation mechanisms.For dense polycrystalline ice, these obstacles are maximal, and the high stresses generated at the crystal boundaries are thought to initiate recrystallization: new crystals are nucleated and crystal boundaries migrate to replace the highly distorted zones (Meyssonnier et al., 2009b;Montagnat et al., 2015).With increasing porosity, the obstacles to dislocation at grain boundaries are successively reduced, and grains can deform more independently from each other with the free space permitted by the pores (Theile et al., 2011).For firn and snow with a density lower than 550 kg m 3 , an additional mechanism is thought to come into play: grain boundary sliding (Alley, 1987).It provides a convenient formalism for modeling firn densification (Lundin et al., 2017;Schultz et al., 2022) but lacks experimental evidence so far (Meyssonnier et al., 2009b).
For snow with even lower density, less attention was drawn to the crystalline nature of ice in snow (e.g., Shapiro et al., 1997).Only recently Riche et al. (2013); Calonne et al. (2017); Montagnat et al. (2021) measured the crystallographic orientations statistically in snow.Indeed, the stress concentration in a few bearing force chains and bonds, related to the shape and size of the ice skeleton (geometrical microstructure), is thought to be the central control of the overall mechanical behavior (Johnson & Hopkins, 2005;Kry, 1975;Wautier et al., 2017).However, the geometrical microstructure alone appears insufficient to describe the diverse creep rates observed on different snow types (Calonne et al., 2020;Fourteau et al., 2022).Besides, the driving deformation mechanism at the microscale, either inter-crystalline deformation (e.g., grain boundary sliding) or intra-crystalline deformation, remains a matter of debate because of a lack of experimental evidence (Meyssonnier et al., 2009a;Sundu et al., 2024;Theile et al., 2011).
An alternative approach to understand and predict the mechanisms at play in snow viscoplasticity relies on numerical experiments.The idea is to evaluate how different hypothetical mechanisms in ice would affect the overall mechanical behavior of snow.Experimental data, even if they do not exhibit the obvious signature of the involved mechanisms, can then be used to discriminate the modeling assumptions.The main ingredients are the porous and crystalline microstructure of snow and the deformation mechanisms active in ice crystals and at their boundaries.Different studies followed this modeling path in the last decades.Johnson and Hopkins (2005); Kabore et al. (2021) modeled the snow microstructure as a set of solid discrete elements interacting with each other through an elastic viscous-plastic contact law.By definition, this method describes the dominant mechanism as grain boundary sliding.The snow viscosity simulated by Johnson and Hopkins (2005); Peters et al. (2021) required an important scaling to reproduce experimental creep data of natural snow.However, the artificial microstructure used in the simulation limits any 1:1 comparison to experimental data.In contrast, Theile et al. (2011); Wautier et al. (2017) used a microstructure directly derived from tomographic data.They only considered the deformation within the crystals and obtained a fair agreement with experimental data.However, Theile et al. (2011) described the snow microstructure as a collection of inter-connected finite element beams, which cannot reproduce deformation obstacles at the crystal boundaries.Wautier et al. (2017) described ice in snow as a homogeneous material, polycrystalline ice.This approach is not suited to account for the free space around ice crystals in snow, which may favor basal dislocations glide, in contrast to dense polycrystalline ice.
In line with these previous studies, we aim to quantify the effect of ice mechanics on snow viscoplasticity.We use time-series of three-dimensional images obtained via tomography during mechanical tests as direct inputs of a mechanical model.The material law for ice in the snow is either homogeneous polycrystalline ice (ice foam model), or each crystal is associated with a given crystal plasticity model (sintered crystal model).This general approach avoids any over-simplification of the snow microstructure and the ice mechanics.

Experimental Data
We used two different mechanical tests: a strain-rate-controlled experiment (Bernard, 2023) and a load-controlled experiment (Bernard et al., 2022) (Figure 1a).In these tests, the macroscopic stress or strain was measured, and the microstructure was regularly scanned at a micro-metric resolution in an X-ray tomograph (DeskTom130, RXSolutions).The experiments were conducted under isothermal conditions on natural snow that was let evolve in a cold room at respectively 20 and 6°C for several weeks.Experimental conditions and imaging setup are summarized in Table 1, and detailed in Bernard et al. (2022); Bernard (2023).
The greyscale attenuation images measured by tomography were binary segmented into pore space and a continuous ice matrix (Bernard et al., 2022;Hagenmuller et al., 2013) (Figure 1b top).This data does not contain any information about crystal boundaries.The ice matrix was thus segmented into individual grains based on two geometrical criteria: negative minimal principal curvature and low contiguity between the grains (Hagenmuller,  Note.The type of snow is defined according to (Fierz et al., 2009), where DF stands for Decomposing and Fragmented snow, and RG stands for Rounded Grains.The numbers shown in brackets, for example, [a, b], mean that the parameter ranges between, for example, a and b.Details can be found in the associated references.
Chambon , et al., 2014;Peinke et al., 2020) (Figure 1b bottom).We assumed that the grains detected by the algorithm correspond to the ice crystals, which is reasonable for this type of snow according to Arnaud, Gay, et al. (1998).

Constitutive Law for Ice
A key ingredient of the numerical model is the elasto-viscoplastic law used for ice.The strain tensor in ice can be decomposed as the sum of an elastic and a viscoplastic part: ɛ = ɛ e + ɛ vp .The elastic strain ɛ e is related to the stress tensor σ and the stiffness tensor C as σ = C : ɛ e .The viscoplastic strain rate εvp is related to the stress tensor σ as εvp = f (σ).The values of C or f depend on whether ice is considered a homogeneous material or a set of sintered single-crystals.

Sintered Crystal Model
Ice in snow can also be explicitly modeled as a set of individual ice crystals sintered together, each characterized by a large viscoplastic anisotropy.Since the crystalline orientation is not captured by tomography, it was randomly sampled in an isotropic distribution (see Text S3 and Figure S4 in Supporting Information S1 for sensitivity analysis).The elasticity tensor C was supposed to be isotropic transverse in the crystal reference frame and defined according to Gammon et al. (1983) (C 11 = 13.9GPa, C 33 = 15.0GPa, C 44 = 3.0 GPa, C 12 = 7.1 GPa, C 13 = 5.8 GPa).For the viscoplastic part, we used the crystal plasticity model of Lebensohn et al. (2009).In this model, ice crystals can deform through slip on three soft basal systems, three hard prismatic systems, and six hard pyramidal systems (Montagnat et al., 2014).On the kth system, the slip-rate γ(k) is related to the shear stress τ (k) by with n (k) the stress exponent, γ(k) 0 the reference shear-rate and τ (k) 0 the critical shear stress.The critical stresses τ (k) 0 for the prismatic and pyramidal systems were assumed to be constant (no hardening) and equal to 260 MPa, which is 20 times larger than the values chosen for the basal systems (13 MPa).For all systems, the creep exponent is n (k) = 3 and the reference shear-rate is γ(k) 0 = 1 s 1 .This setup enables to reproduce the behavior of polycrystalline ice, when the single crystals are randomly arranged into a dense packing (Text S1 in Supporting Information S1).This model, although simpler than Suquet et al. (2012), was preferred for the sintered crystal model, as it allows a simple comparison with the foam model, giving the same stress exponent n on the macroscopic scale.
To account for temperature effects on ice viscoplasticity, we scaled the pre-factor A of the ice foam model (Equation 1) and the pre-factor γ(k) 0 of the crystal model (Equation 2) by an Arrhenius relation: e Q RT , with Q = 69.1 kJ mol 1 the activation energy of ice, R = 8.3 J (mol K) 1 the universal gas constant and T the temperature (K) (Mellor & Testa, 1969).This scaling enables us to account for the different temperatures in the two experiments ( 8 and 18.5°C, Table 1).

Simulation Set-Up
The simulations were performed with the Fast Fourier Transform (FFT)-based solver AMITEX_FFTP (Gélébart et al., 2020).This solver takes as direct input the 3D air-ice image (ice foam model) or the 3D image of the crystal assembly (sintered crystal model).The numerical model cannot account for topological change of the ice skeleton (e.g., new contacts).The simulation of one experiment thus consists of several simulations of "instantaneous creep experiments" as in Theile et al. (2011), with initial conditions varying through the test (re-set at each new image).The Norton-Hoff model and the crystal plasticity model were implemented through the MFront code generator (Helfer et al., 2015).The numerical integration of crystal plasticity is computationally expensive, but the solver benefits from the MPI implementation.
Using an FFT-based solver involves using periodic boundary conditions.The boundary conditions were chosen to mimic the mechanical tests, where the snow sample cannot deform laterally and the axial friction on the sample holder is limited.Therefore, in the simulation, we imposed the average lateral deformation (E xx = E yy = 0) and the shear stress to zero (Σ xy = Σ yz = Σ xz = 0), where E and Σ are respectively the volume average of strain and stress.The axial loading was imposed according to the type of test (Figure 1c).
To reduce the computation time, the simulations were not conducted on the entire scanned volume at the nominal image resolution.Instead, we reduced the resolution by a factor of 2, that is, to voxels of size 16 and 17 microns respectively for the strain-rate-controlled and load-controlled experiments, and we extracted a volume of 300 × 300 × 300 voxels inside the scans (Figure 1).We performed a volume convergence analysis on one sample and showed that the chosen volume and resolution are representative of the sample.To fasten the solver convergence, we reduced the infinite mechanical contrast between ice and air to E ice /E air = 10 8 by giving a very small stiffness to air.We observed that this choice reduced the computing cost and did not affect the results (Figure S2 in Supporting Information S1).

Simulated Macroscopic Stress-Strain Curve
We observe three different regimes of stress evolution with strain during the compression test at a constant strainrate (Figure 2a inset).First, the stress increases linearly with strain in an elastic regime characterized by a slope E.Then, the stress deviates from elasticity with viscoplasticity progressively activated.After this transient regime, the stress eventually reaches a constant value, the yield stress σ Y , where the snow sample perfectly flows (stationary creep).This yield stress depends on the constitutive law of ice and the snow microstructure.Similarly, Figure 2b (inset) shows the time evolution of strain under constant stress.After a brief phase of pure elasticity during the numerical loading of the sample (stress from 0 to 2.1 kPa in 1 s), the ice matrix starts to flow and the strain rate increases to a constant value.
The permanent viscoplastic regime is reached after a typical strain of 0.1% (Figure 2a) or a time of 10 4 s (Figure 2b).The yield stress or flow rate are, hereafter, derived at a simulation strain of 0.5% and a simulation time of 10 5 s.The time or strain required to reach the stationary creep flow regime is relatively small compared to the time between two scans (∼8 hr) or change in the microstructure.Therefore, the assumption that mechanical creep experiments with an evolving microstructure can be decomposed into "instantaneous creep experiments" seems reasonable.In other words, this separation of time scales enables us to compare the yield stress estimated numerically on a given microstructure to the stress measured when the microstructure was captured (Theile et al., 2011).
The elastic and viscoplastic regimes were simulated with different models for ice.The stiffness of the sample is not affected by the modeling assumptions, that is, the slight anisotropy of crystalline elasticity (30%) could be neglected.Indeed, a relative error of 0.5% is committed on the axial Young's modulus with the ice foam model compared to the sintered crystal model.Therefore, assuming ice in snow as a foam of ice to study snow elasticity (Hagenmuller, Calonne, et al., 2014;Hagenmuller, Theile, & Schneebeli, 2014;Köchle & Schneebeli, 2014;Reuter et al., 2019;Schneebeli, 2004;Srivastava et al., 2010;Wautier et al., 2015) appears reasonable, even never properly evaluated so far.In contrast, the simulated viscoplastic regime is very different between the ice foam and the sintered crystal models.For instance, at a strain rate of 1.8 × 10 6 s 1 , a yield stress of 22.2 and 8.7 kPa is obtained respectively for the ice foam and sintered crystal models (Figure 2a).As expected, the crystal model makes the snow sample flow more easily.Interestingly, blocking the hard slip systems of the sintered crystal model (only basal glide) does not significantly affect the macroscopic mechanical behavior of the snow for both experiments.
As the snow stress strain-rate relation is non-linear, the classical compactive viscosity cannot describe both tests (Kojima, 1967).The creep exponent of all slip systems is here constant (n = 3).Therefore, per homogenization, this exponent scales to the visco-plastic model of snow (Wautier et al., 2017) and we can identify the observed creep of snow by a Glen's law with a non-linear viscosity B:

Evaluation of the Numerical Model With the Experiments
In the strain-rate-controlled experiment, the measured non-linear viscosity B increases from 6.8 × 10 2 to 1.2 × 10 4 MPa n s 1 , while the density increases from 206 to 458 kg m 3 (Figure 2).The simulations well reproduce this trend with density, and the order of magnitude of B. For a similar density range, Scapozza and Bartelt (2003) obtained a viscosity change from 10 2 -10 6 MPa n s 1 at a temperature of 12°C.Note that in this experiment, the exponent of Glen's law was identified experimentally, and varies from 1.7 for a density of 150 kg m 3 to 3.7 for a density of 423 kg m 3 .In the load-controlled experiment, the increase of the non-linear .The inset figures correspond to the simulated macroscopic stress or strain evolution for different constitutive materials for a given microstructure.The error bars are related to the limited simulation volume (error of ±14.9%, Figure S3 in Supporting Information S1), the resolution decrease (error ±16.4%, Figure S3 in Supporting Information S1) and uncertainties on temperature (error of ±6.2% at 8 and ± 6.8% at 18.5°C).

Geophysical Research Letters
10.1029/2023GL107676 viscosity with time is less pronounced, as the overall deformation or densification remains rather small (ϵ z ≤ 9.1%).However, the simulation also reproduces this tedious trend driven by microstructure evolution.
The simulations with the sintered crystal and the ice foam model respectively underestimate and overestimate B in the strain-rate-controlled experiment.The measured and simulated value of B at the beginning of the experiment is uncertain due to the limited resolution of the force sensor (±4 kPa).However, for stresses higher than 25 kPa, we observe that the sintered crystals model is closer to the experimental measurements, than the foam model.Indeed, the mean absolute percentage error (MAPE) between simulations and the experiments is equal to 414% for the ice foam model and to 62% for the sintered crystal model.The difference between the sintered crystal and the ice foam models reduces when density increases.The ratio decreases from 16.6 at a density of 206 kg m 3 , to 9.7 at a density of 458 kg m 3 .This is consistent with the fact that the two models converge together when the porous structure becomes pure ice (see Supp.Inf.S1 in Supporting Information S1 and R. Lebensohn (2001)).
More generally, we expect that the models become closer when the crystal boundaries become pre-dominant barriers to basal dislocation glide.
For the load-controlled experiment, the simulations based on the sintered crystal model are closer to the experiment (over-estimation of B by an average factor of 8.9) compared to the ones based on the ice foam model (over-estimation of B by 157).Note that the error factor estimated here on the non-linear viscosity are raised to the power of 3 compared to the errors computed on linear viscosity (η = σ ε ) , as it is presented by (Kojima, 1967;Theile et al., 2011).Experimental uncertainties also exist in the load-controlled test (Figure 2b), but they are difficult to estimate and related to (a) the measurement of the strain rate and (b) the force effectively applied to the considered snow volume.Nevertheless, they can be estimated to be less than the ratio 8.9 between the best model and the experiment.
Overall, for snow, here with a density in the range 206-458 kg m 3 , the two models for ice yield very different viscoplastic properties, and the sintered crystal model appears to be closer to the experimental data.This observation agrees with the recent work of Fourteau et al. (2023), who shows that it is impossible to find a calibrated isotropic law to represent different types of microstructures, and questions the magnitude of non-linear viscosity simulated by Wautier et al. (2017).Nevertheless, the significant discrepancies between the sintered crystal model and experiments show that some mechanisms may still be missing in this modeling.

Microscale Drivers of Snow Viscoplasticity
In the sintered crystal model, we observe that a small deformation of the ice matrix yields important deformation of the pore space and the overall snow material.Indeed, the average deformation of the ice matrix is here typically 100 times smaller than in the equivalent snow material (Figure S4 in Supporting Information S1).In addition, the distribution of strain or stresses in the ice matrix itself is highly heterogeneous (standard deviation of 30% on strain in the ice matrix) and comprises both zones in tension and compression (noted with a positive sign) (Figure 3).For instance, 4% of the ice volume exhibits an axial deformation 10 times larger than the average ice deformation, while a significant part (e.g., 34.5%) of the microstructure exhibits almost no deformation (e.g., |ϵ z | ≤ 0.1ϵ z ).This confirms previous numerical observations on elastic properties (Hagenmuller, Theile, & Schneebeli, 2014) but is here even more pronounced as the ice can locally undergo large viscoplastic deformation.
Interestingly, similarly to dense ice, the basal slip system contributes at most to the total simulated plastic deformation of snow (90% in average for ice, 99% for snow, see Figure S4 in Supporting Information S1).However, the entanglement of slip system deformations appears to be different between snow and dense ice.The model shows that the deformation accommodation between crystals is permitted by the pore space and does not require prismatic and pyramidal slips, whereas the latters are pre-requisite for the simulation of any deformation in dense ice.In addition, blocking the hard slip systems of the sintered crystal model does not significantly affect the macroscopic mechanical behavior of the snow for both experiments (Figure 2).Thus, a stronger effect of the basal system parameterization on the macroscopic response of the snow can be expected.Indeed, the chosen parameterization was mostly evaluated on pure ice experiments (Castelnau et al., 1996(Castelnau et al., , 1997;;Mansuy et al., 2000), but different choice for the basal slip system (e.g., exponent n or critical shear stress τ 0 ) might yield similar results for polycrystalline ice but not for snow (Duval et al., 1983a;Suquet et al., 2012).
Stresses and strains are heterogeneously distributed in the microstructure and generally located at necks between ice grains that participate in a force bearing chain (Figure 3).These geometric necks also mainly correspond, here, to the crystal boundaries.The grain boundaries (defined as a one voxel wide band on either side of the interface between 2 grains, which represent 11% of the ice volume fraction), concentrate 20.5% of the equivalent stress, 21% of basal equivalent strain and 22.2% of non-basal equivalent strain.When we do not account anymore for the single crystals (i.e., foam model), the stress concentration reduces to 17% of the equivalent stress.This stress concentration is thus not only linked to the geometry of the microstructure, but also to the mechanical model of ice.Interestingly, the relative basal activity remains high everywhere.It is lower next to certain crystal boundaries, where the stresses are the highest (Figure 3d).Non-basal contributions also appear in crystals with very little plastic deformation but large stress concentration (e.g., the grain circled in Figure 3).This observation indicates that crystal boundaries may slightly limit basal slip, but in a proportion that is not visible on the overall viscous behavior (here for a snow sample with a density of 235 kg m 3 ).

Conclusion
We simulated the viscoplasticity of snow based on its 3D microstructure and different constitutive laws for the ice matrix.These numerical experiments were compared to cold-room in-tomograph experiments, either loadcontrolled or strain-rate-controlled.The macro-scale comparison revealed that ice in snow rather behaves as a set of sintered crystals than a foam of polycrystalline ice.Moreover, we showed that the accommodation of deformation by means of the hard non-basal slip systems is hardly needed, even less than for bulk polycrystalline ice.This is attributed to the ability to relax strain incompatibilities at ice/air interfaces.The residual mismatch between the measured and the simulated viscosity tends to demonstrate that other mechanisms occurring, for example, at bonds need to be accounted for, such as, role of non-basal contributions with hardening (Duval et al., 1983a;Suquet et al., 2012), superplasticity (Alley, 1987;Goldsby & Kohlstedt, 2001;Raj & Ashby, 1971;Sundu et al., 2024), but also ductile failure (Capelli et al., 2020;Kirchner et al., 2001).

Figure 1 .
Figure 1.Main workflow.(a) Schemes of the mechanical tests.(b) 3D snow microstructure image used for the simulation: binary ice-air image for the ice foam model (top), grain image for the sintered crystal model (bottom).(c) Numerical set-up and boundary conditions.

Figure 2 .
Figure2.Comparison of the experimental and numerical results.(a) Evolution of the nonlinear viscosity and the yield stress as a function of strain (strain-ratecontrolled).(b) Temporal evolution of the nonlinear viscosity and the strain rate (load-controlled)).The inset figures correspond to the simulated macroscopic stress or strain evolution for different constitutive materials for a given microstructure.The error bars are related to the limited simulation volume (error of ±14.9%,FigureS3in Supporting Information S1), the resolution decrease (error ±16.4%, FigureS3in Supporting Information S1) and uncertainties on temperature (error of ±6.2% at 8 and ± 6.8% at 18.5°C).

Figure 3 .
Figure 3. Slice view of the tested sample shown in Figure 1 and associated mechanical fields: (a) Von Mises equivalent stress σ eq .(b) Equivalent viscoplastic basal slip.(c) Equivalent viscoplastic non basal slip.(d) Basal activity.Strain-rate-controlled experiment, axial stress corresponding to this microstructure of 8.7 kPa.

Table 1
Experimental Conditions of the Mechanical Tests