Numerical homogenization of the viscoplastic behavior of snow based on X-ray tomography images

While the homogenization of snow elastic properties has been widely reported in the literature, homogeneous rate-dependent behavior responsible for the densification of the snowpack has hardly ever been upscaled from snow microstructure. We therefore adapt homogenization techniques developed within the framework of elasticity to the study of snow viscoplastic behavior. Based on the definition of kinematically uniform boundary conditions, homogenization problems are applied to 3-D images obtained from X-ray tomography, and the mechanical response of snow samples is explored for several densities. We propose an original postprocessing approach in terms of viscous dissipated powers in order to formulate snow macroscopic behavior. Then, we show that Abouaf models are able to capture snow viscoplastic behavior and we formulate a homogenized constitutive equation based on a density parametrization. Eventually, we demonstrate the ability of the proposed models to account for the macroscopic mechanical response of snow for classical laboratory tests.


Introduction
It is now well known that the macroscopic mechanical behavior of snow strongly depends on its microstructure, which is mainly characterized by its density (Mellor, 1974) and its topology (Shapiro et al., 1997), by the mechanical behavior of the ice matrix 10 (elastic, visco-plastic, brittle-failure) depending on the external load, by the temperature (Schweizer and Camponovo, 2002) and the applied strain rate (Schulson et al., 2009). For example, at high strain rate, large deformations of snow are mainly controlled by grain rearrangements resulting from the failure of cohesive bonds, whereas at very low strain rates (< 10 −5 s −1 typically), snow exhibits a visco-plastic behavior (Narita, 1984;Salm, 1982) which plays an important role in the long-term densification of the snowpack. In practice, a good knowledge of the macroscopic mechanical behavior of snow in a wide range 15 of applied loads, strain rates and temperatures is of a particular interest with respect to avalanche risk forecasting or to determine the forces on avalanche defense structures.
During the last decades, numerous experimental studies have been performed in order to characterize the macroscopic behavior of different types of snow under various loading conditions and temperatures (Mellor, 1974;Salm, 1982;Desrues et al., 20 1980;Shapiro et al., 1997;Bartelt and von Moos, 2000;Moos et al., 2003;Scopozza and Bartelt, 2003a). In the framework of the continuum mechanics, several models have been then proposed in order to reflect these experimental data (Desrues et al., 1980;Scopozza and Bartelt, 2003b;Cresseri and Jommi, 2005;Navarre et al., 2007;Cresseri et al., 2009). However, the fitted material parameters arising in these models often characterize the mean properties of a few types of snow in a restricted density range. Thanks to the recent application of X-ray tomography to snow (Brzoska et al., 1999;Schneebeli, 2004;Kaempfer et al., 2005;Flin and Brzoska, 2008;Chen and Baker, 2010;Srivastava et al., 2010;Pinzer et al., 2012;Wang and Baker, 2013;Adams and Walters, 2014;Calonne et al., 2015) good databases of 3D images for the different snow types described 5 in the international classification (Fierz et al., 2009) are now available (Calonne et al., 2012;Löwe et al., 2013). Given these extensive geometrical descriptions of snow, its corresponding macroscopic behavior can be up-scaled in a more systematic way thanks to the use of techniques derived from the homogenization theory (Dormieux and Bourgeois, 2002;Auriault et al., 2010).
In recent works, the combination of X-ray tomography imaging, finite elements techniques or discrete element methods (DEM) 10 and ever increasing computing power was used to bridge the gap between the topology of the ice skeleton of snow and its mechanical behavior. If discrete element methods are well suited to model deformations resulting from grain rearrangement and grain bound breakage (Hagenmuller et al., 2015), finite element modeling is more adapted to account for the deformation of the ice skeleton. This last method is thus well-suited to the long term compaction of the snowpack under its own weight coming from the ice viscosity. However, most of the studies only considered the elastic behavior of snow so far (Schneebeli, 2004;15 Pieritz et al., 2004;Srivastava et al., 2010;Köchle and Schneebeli, 2014;Wautier et al., 2015;Srivastava et al., 2016), possibly up to a brittle failure (Hagenmuller et al., 2014). Concerning the modeling of more complex snow constitutive behaviors, the proposed approaches mainly focus on the modeling of uniaxial compression tests. For instance, Theile et al. (2011) proposed a beam network model based on 3D images to simulate creep of snow, Chandel et al. (2014) used an elasto-plastic constitutive law for ice, Schleef et al. (2014) considered the viscous behavior of snow in the case of uniaxial compression tests followed by 20 a relaxation phase.
In the present work, we propose to determine the 3D macroscopic viscoplastic constitutive law of snow in the wake of the homogenization approach performed in Wautier et al. (2015) within the framework of elasticity. For that purpose, typical kinematically uniform boundary conditions (KUBC) are applied to 3D images of snow, and finite element simulations are run 25 in order to link the macroscopic stress response of different snow samples to imposed strain rates. The viscous mechanical behavior of ice is described by a power law of exponent n as in Theile et al. (2011). Due to the non-linear behavior under consideration, the homogenization does not provide the complete structure of the macroscopic constitutive equation (Auriault et al., 1992(Auriault et al., , 2002Geindreau and Auriault, 1999;Orgéas et al., 2007). However, it can be shown that the exponent n is preserved at the macroscopic scale and the macroscopic dissipation potential is the volume-averaged of the local dissipation one 30 (Suquet, 1993). Using these properties, the macroscopic constitutive equation of snow can be formulated within the framework defined by the theory of tensor function representation and by using macroscopic isodissipation surfaces (Green, 1972;Abouaf, 1985;Geindreau et al., 1999b;Orgéas et al., 2007).
The paper is organized as follow. In section 2, the numerical homogenization procedure used in Wautier et al. (2015) is recalled and adapted to the study of a non-linear constitutive equation. The section 3 presents the post-processing procedure that was used in order to characterize the macroscopic viscous behavior of snow in terms of macroscopic isodissipation curves. In section 4, the obtained numerical results for snow samples of different densities are presented and the ability of an Abouaf's model (Abouaf, 1985) to describe the macroscopic viscous behavior of snow is highlighted. Finally, in section 5, the mechan-5 ical responses of snow for classical experimental tests (uniaxial, oedometric and triaxial compression tests) are analysed and compared using the proposed model.
2 Numerical homogenization procedure: from image to macroscopic mechanical response Based on the homogenization theory, it is often possible to replace a heterogeneous material by an equivalent homogeneous one provided that its microstructure is sufficiently small with respect to the macroscopic scale of interest. With respect to snow, 10 this separation of scale hypothesis is satisfied in most of the cases and its macroscopic mechanical behavior can be deduced from mesovolumes obtained thanks to X-ray tomography. Previous studies showed that in most of the cases, samples of a few centimeters can be considered as representative elementary volumes (REV) for the study of the mechanical behavior of snow (Wautier et al., 2015;Srivastava et al., 2016). 15 Irrespective of the size of the sample considered, the boundary conditions used in an homogenization procedure introduces undesired boundary effects of varying thickness. Depending on the type of boundary conditions used, the size of the REV should be adapted accordingly. Three particular types of boundary conditions are considered to give relatively small REV. In decreasing order of REV (Kanit et al., 2003), these are statically uniform boundary conditions (SUBC), with a macroscopic homogeneous stress imposed on the boundary, kinematically uniform boundary conditions (KUBC), with a macroscopic ho-20 mogeneous strain imposed on the boundary, and periodic boundary conditions (PBC), with a periodicity condition imposed on the displacement field and the normal stress across the sample boundaries. Although PBC are considered to give the best convergence with respect to the size of the REV (Kanit et al., 2003), their application to a non-periodic highly porous microstructure is not straightforward. It is necessary, for example, to enclose the sample by a virtual boundary or assume that the pores are filled by a soft material. In order to avoid the introduction of such artefacts, KUBC were retained as in Wautier et al. 25 (2015). The KUBC numerical homogenization procedure introduced in Wautier et al. (2015) is used in this paper and easily adapted to the study of the elasto-viscoplastic behavior of snow. It consists in the four steps recalled in Fig. 1.
The first two steps remain unchanged and consist in: (i) meshing of 3D-microtomographic images (Step 1), (ii) defining the kinematic relation u = E · x between the homogeneous macroscopic strain E and the displacement field u on the boundary 30 ( Step 2).
Step 1 requires the use of the MATLAB open-source toolbox iso2mesh (Fang and Boas, 2009) while step 2 is achieved thanks to the use of the plug-in Homtools (Lejeunes et al., 2011). More details can be found in Wautier et al. (2015). The next two steps are modified in order to take into account the change in the constitutive modeling of ice.

Abaqus plugin HomTools
Step 4 Step 3 Figure 1. Four-step procedure used in order to transform 3-D microtomograph images of snow into finite element models and numerically solve KUBC homogenization boundary value problems.

Elasto-viscoplastic behavior of ice (Step 3)
In the following, the mechanical behavior of the polycrystalline ice is supposed to be elasto-viscoplastic and isotropic. The total strain rate tensor (ε) is decomposed as the sum of the elastic part (ε e ) and a viscous part (ε v ) aṡ The elastic part can be expressed as, ε e = (C ice ) −1 : σ, where C ice is the elastic stiffness tensor, σ is the Cauchy stress tensor 5 and ":" the double contraction product. Due to isotropy, C ice is fully defined by a Young's modulus E and a Poisson ratio ν.
Concerning the viscous part, at low strain rates, the non linear mechanical behavior of ice is usually described by a power law (Mellor, 1974;Schulson et al., 2009), i.e. the Norton Hoff law in 3D (Lemaitre and Chaboche, 1985). Thus, we havė where A and n are two material parameters (which usually depend on the temperature), σ is the deviatoric stress tensor and σ eq (σ) is the equivalent stress defined as where I is the second order identity tensor, Tr is the trace operator ands 2 is the second invariant of σ. It can be shown (Lemaitre and Chaboche, 1985) that the viscous strain rate tensorε v derives from a viscous potential ω(σ) as From (2) and (4), one can define the equivalent strain rateε eq (ε v ) as the dual variable of the equivalent stress σ eq such thaṫ whereē 2 is the second invariant of the deviatoric part ofε v .
If p v stands for the volumetric mechanical dissipation, the equivalent stress and the equivalent strain rate verify At the microscopic scale, the viscoplastic deformation of ice is incompressible. Consequently the equivalent stress (resp. the equivalent strain rate) depends only on the second invariants 2 of σ (resp. the second invariantē 2 ofε v ).
Overall, the ice matrix is thus modeled as an isotropic elasto-visco-plastic material in the finite element commercial software Abaqus. Let us remark, that a pure visco-plastic behavior cannot be simulated in a Lagrangian code such as Abaqus and the 15 knowledge of the ice Young's modulus (E) and Poisson's ratio (ν) are also needed to ensure the efficiency of the numerical scheme. The values of the four material constants in this constitutive modeling, namely A, n, E and ν, were determined thanks to a uniaxial compression test on a cylinder of polycrystalline ice. The experiment was carried out at -10 • C on an ice cylinder thanks to a micro-press operating in cold-room. After a first compression step realized at a constant strain rate of 2 × 10 −6 s −1 until reaching a stress stabilization around 475 kPa, the ice viscous properties were identified thanks to a relaxation test. After 20 3500 s and as the axial stress was reduced by half, the ice elastic properties were deduced from the unloading of the ice cylinder at a higher strain rate of 1 × 10 −3 s −1 . The Table 1 summarizes the parameter values used in this paper, which are consistent with the ones in the literature (Chandel et al., 2014;Schulson et al., 2009).

Macro-strain paths definition (Step 4)
Given a time dependent macroscopic strain loading E(t), the KUBC homogenization problem to be solved reads as where V stands for the domain occupied by the whole snow sample and ∂V its boundary. The spatial heterogeneity of the mechanical properties of snow is captured thanks to the functions A(x) and E(x) defined as where V i ⊂ V is the domain occupied by the ice matrix. Similarly to the elastic case (Wautier et al., 2015), the macroscopic stress tensor Σ is deduced from the knowledge of its microscopic counterpart thanks to the volume averaging As a result, for a given macroscopic strain loading E(t), the macroscopic stress response Σ(t) is recovered. The implicit 10 function linking these two second order tensors characterizes the homogeneous behavior of the snow sample considered and can be put in the form: whereĖ e is the macroscopic elastic strain rate tensor andĖ v is the macroscopic viscous strain rate tensor. The elastic part can be expressed as, E e = (C hom ) −1 : Σ, where C hom is the homogenized stiffness tensor (Wautier et al., 2015). This tensor can be obtained by performing only six simulations on Representative Elementary Volumes extracted from 3D images. By contrast, the homogenization of the visco-plastic behavior requires a priori an infinite number of numerical simulations. However, this number of simulations can be reduced by taking into account some theoretical results (Auriault et al., 1992;Suquet, 1993;5 Orgéas et al., 2007). Indeed, it can be shown that: -The homogeneity of degree n of the microscopic viscous constitutive equation (2) is preserved in the homogenization process. In other words, the macroscopic viscous strain rateĖ v is a homogeneous function of degree n of the macroscopic stress Σ, and the macroscopic volumetric mechanical dissipation P v =Ė v : Σ is an homogeneous function of degree n + 1 of Σ: As a result, the choice in the macroscopic strain rateĖ v can be reduced to the unit sphere in the second order tensor space.
-The macroscopic dissipation potential Ω(Σ) is the volume-average of the local dissipation potential ω 15 and consequently, as at the microscopic scale (see equation (4)), we have: where Σ eq (Σ) is the macroscopic equivalent stress, which verifies withĖ eq (Ė v ) the macroscopic equivalent strain rate defined by duality. As a result, the macroscopic viscoplastic law 20 (13) is perfectly defined if the macroscopic equivalent stress Σ eq is known. The latter equation (14) shows that this macroscopic equivalent stress Σ eq can be fitted on iso-mechanical dissipation surfaces in the space associated with Σ.
In the case of general anisotropy, the form of Σ eq can be formulated within the framework defined by the theory of representation of anisotropic tensor functions (Smith, 1971;Liu, 1982). It is also important to mention that for the ice matrix, the overall response of snow is insensitive to the sign of Σ as a consequence of definition (4). This condition may 25 be expressed as Ω(Σ) = Ω(−Σ). Finally, let us remark that by definition (see (13) and (14)), the macroscopic strain ratė E v is normal to iso-mechanical dissipation surfaces (normality rule). In the following, for the sake of simplicity, we will suppose that the macroscopic viscoplastic behavior of snow is isotropic. In this particular case, it can be shown (Abouaf, 1985;Geindreau et al., 1999b;Danas et al., 2008) that the macroscopic equivalent stress is written: where φ is the snow porosity and (S 1 ,S 2 ,S 3 ) are the three invariants of the macroscopic stress tensor Σ defined as: Similarly, the macroscopic equivalent strain rateĖ eq takes the form: where (E 1 ,Ē 2 ,Ē 3 ) are the invariants of the strain rate tensorĖ v defined as: In contrast with the microscopic scale (see equation (3)), the macroscopic equivalent stress (15) depends on the three invariants (S 1 ,S 2 ,S 3 ) of the macroscopic stress tensor Σ. Indeed, at the macroscopic scale, the viscoplastic snow deformation is compressible. This compressibility, characterized by E 1 , depends on the level of the mean pressure (S 1 /3) applied on the snow sample, as well as the mean shear stress (S 2 ). The third invariantS 3 characterizes the loading type and is linked to the Lode angle θ in the stress space (Lemaitre and Chaboche, 1985;Danas et al., 2008): As a first approximation, it seems reasonable to assume that the influence of the third invariantS 3 is negligible (Green, 1972;Abouaf, 1985;Geindreau et al., 1999b;Fritzen et al., 2012). Consequently, the macroscopic volumetric mechanical dissipation P v depends on the first and second stress and strain invariants and not only on the second ones as at the microscale (6). 20 The relation (20) shows that, for a given snow porosity, the equivalent macroscopic stress Σ eq can be fitted on iso-volumetric mechanical dissipation curves in the plane (S 1 /3,S 2 ). These isodissipation curves can be obtained by plotting the values (S 1 /3,S 2 ) corresponding to different loading conditions defined by (E 1 ,Ē 2 ). Therefore, the choice was made to run numerical simulation for seven diagonal strain rate tensors defined such that the loading direction in the plane (E 1 ,Ē 2 ) varies from 0 • to 90 • . More explicitly,Ė applied on the sample is taken as Finally, to be consistent with the isotropy hypothesis, numerical simulations have been performed on the most isotropic snow samples with respect to their elastic behavior from the snow database used in Wautier et al. (2015). With reference to the supporting information of the cited paper (Wautier et al., 2015), the name and the principal characteristics of each sample are 5 recalled in Table 2. The porosities of the selected samples vary from 0.43 to 0.87, which covers almost the entire range of porosity of seasonal snow. Each sample presents similar correlation lengths ( 1 , 2 , 3 ) (Löwe et al., 2013;Calonne et al., 2014) in the three space directions. All the simulations have been performed on volumes extracted from the 3D images sufficiently large to be considered as REV, as in Wautier et al. (2015). 3 Post-processing procedure: from macroscopic stress response to a homogenized model for snow visco-plasticity 10 From the homogenization procedure presented in the previous section, the time response of a given isotropic snow sample is recovered for the seven loading directions in the plane of the strain invariants (E 1 ,Ē 2 ) given by the equation (21). In order to deduce the overall viscous behavior of the snow considered from these few tests, a post-processing procedure consisting in basically three steps is proposed and summarized in Figure 2.
3.1 Extracting the viscous response (step a) 15 Because snow is locally modeled as an elasto-visco-plastic material in Abaqus (see subsection 2.1), the macroscopic time response Σ(t) = 1 |V | V σ(x) dV deduced from numerical simulations is composed of an elastic and a viscous part. The typical mechanical response of a snow sample to a constant strain rate is illustrated on Figure 3 for the snow sample RG_1600 (Table 2) under a strain rate characterized by a loading angle θ = 65 • (equation (21)). In order to reckon the influence of the elastic constitutive behavior on the overall mechanical response, the time response of the two first stress invariants S 1 and 20S 2 is plotted for two different Young's moduli. In both cases, the initial mechanical response of the sample is driven by the elastic properties with a sharp increase in S 1 andS 2 proportional to the Young's modulus under consideration (Figure 3).
Step a: -Compute the mechanical response . for a given constant strain rate -Determine the viscous response associated with the strain state Step b: -For each loading given in Eq. (21), compute the associated viscous dissipated power with Eq. (22) -Rescale the numerical points on an isodissipation curve thanks to Eq. (23) Step c: -Fit an Abouaf's model thanks to Eq. (24) and (27) Figure 2. Three-step post-processing procedure used in order to formulate a homogenized viscous constitutive equation.
Then a transient regime is observed followed by a permanent one in which the mechanical response is mainly driven by the viscoplastic behavior of ice. In the case where the ice viscosity is activated everywhere in the ice skeleton, the macroscopic stress Σ(t) should stabilize around a constant value according to equation (2) and to the fact that the strain rate is taken constant on average (see equation (21)). However, both stress invariants S 1 andS 2 follow a non horizontal linear asymptote which slope is proportional to the Young's modulus value. Indeed, when E is increased by a factor 2, the slope of the linear asymptote is 5 multiplied accordingly (Figure 3).
In order to extract the viscous component of the stress in the permanent regime, the initial an final asymptotes are plotted in  Table 2) under a strain loading characterized by θ = 65 • in equation (21). Two Young's moduli are used and the resulting time responses are compared.

Computing isodissipation curves (step b)
For a given snow sample of porosity φ and for each applied loading path (21), the macroscopic volumetric mechanical dissi- (20) is computed as where S v 1 andS v 2 are the characteristic stress invariants obtained in the first step of the post-processing process (Figure 2).

5
Each loading path leads to different values of P v . However, iso-mechanical dissipation points in the plane (S 1 /3,S 2 ) can be recovered thanks to the homogeneity property (11). Given an arbitrary value of P • v = 1 Pa.s −1 , the corresponding macroscopic strain and stress invariants are computed as 11 The Cryosphere Discuss. Thanks to this rescaling, the seven homogenization tests (21) enable the description of an isodissipation curve in the plane of the stress invariants (S 1 /3,S 2 ) as illustrated in Figure 2 (step b). For each point (S • 1 /3,S • 2 ) on this graph, the associated flow vector (E • 1 ,Ē • 2 ) is plotted. The viscous dissipated power is thus simply equivalent to the scalar product S

Abouaf's model (step c)
Within the framework presented in section 2.2, Abouaf (1985) has suggested to use the macroscopic equivalent stress proposed 5 by Green (1972) to describe the viscoplastic behavior of metal powders at high temperatures. This macroscopic equivalent stress Σ eq (Σ) is written where f (φ) and c(φ) are two material functions which depend on snow porosity only. When φ = 0, we have f (φ) = 0 and c(φ) = 1 in order to recover the equivalent viscous stress of the ice matrix (3): Σ eq (Σ, φ = 0) = σ eq (σ). From the definition 10 of the viscous strain in (13) together with the previous definition of the equivalent stress in (24), it can be shown thaṫ As a result, the corresponding macroscopic equivalent strain rate introduced in (20) readṡ For a given porosity φ, the combination of (24) and (20) provides an implicit definition of f (φ) and c(φ) such that, for all In practise, the optimal values for f (φ) and c(φ) are obtained by minimizing the quadratic error between the model and the numerical points.
4 Results and discussion 20 The homogenization and the post-processing procedure presented in the previous sections are applied to six isotropic snow samples of various densities chosen in the same database as Wautier et al. (2015) and already introduced in Table 2. In Figure   4, the seven points (S • 1 /3,S • 2 ) corresponding to the strain rates of equation (21) Table 2. It should be underlined here that each isodissipation curve is typical of a given snow characterized by its density and thus each curve is also an iso-density curve.  Figure 4. Isodissipation curves in the plane of the stress invariants (S1,S2) corresponding to an arbitrary dissipated power P • v for the six snow samples of Table 2

Isodissipation curves for different snow samples
The comparison between the simulated points and the Abouaf's model shows overall good agreement in terms of stress prediction even if a systematic slight deviation is observed between the model and the simulated points for the highest S 1 values.
It should also be noticed that the stress state corresponding to an isotropic strain rate (Ē 2 = 0) is not completely isotropic (S 2 = 0). This feature cannot be captured by the Abouaf modeling, which assumes a perfect isotropic behavior of the material.

5
Even if the snow samples used in this study were selected as isotropic as possible, a slight anisotropy should account for the observed residual deviatoric stress component existing under an isotropic strain loading.
As already mentioned in section 2.2, the viscous behavior of snow should be insensitive to the sign of Σ as the ice matrix behave exactly the same in tension and in compression. In the stress space (S 1 /3,S 2 ), this results in the symmetry of the isodissipation curves with respect to the axis S 1 /3 = 0. Provided that the isodissipation curves are smooth, their tangent for 10 S 1 /3 = 0 is horizontal which is respected in Figure 4. It must be mentioned that when snow is subjected to large strain levels, geometrical effects will introduce non linear effects and the mechanical response in tension will differ from the one in compression. These effects can also be investigated using the same homogenization procedure.

Density dependence of the isodissipation curves
As the snow density increases, the isodissipation curves tend to expand, and conversely the flow vectors tends to get smaller. 15 In terms of physics, this means that the denser the snow, the smaller the applied strain rate in order to dissipate the same level of viscous power. In the meantime the applied stress should be increased. This is consistent with the fact that fresh snow tends to get denser more rapidly than already compacted snow under the same imposed loading.
The density dependence of the snow viscous behavior is fully described by the evolution of f (φ) and c(φ) with respect to the snow compacity (ρ snow /ρ ice = 1 − φ) represented in Figure 5. As the snow density decreases in Figure 5, the parameter values 20 increase, which is consistent with the implicit definition of the isodissipation curves in equation (24). For a given equivalent stress Σ eq , higher values for f and c will result in lower stress invariants S 1 andS 2 as observed in Figure 4.
As classically done for metal powders (Geindreau et al., 1999b), some phenomenological fits can be proposed for f and c.
In order to account for the change in the porosity range between metal powders and snow, the compacity limit value of 0.57 introduced in Table II in Geindreau et al. (1999b) is set equal to zero. As a result, the general form of the proposed fits is written The above fits respect the theoretical values f (0) = 0 and c(0) = 1 already mentioned in section 3.3. For a highly porous snow (φ → 1), an infinitely small stress level would be needed in order to produce a high viscous dissipation. This is consistent with the infinitely high values for f and c proposed by the above phenomenological fits (28). These fits allow a good description 30 of the numerical points resulting from the homogenization of the six snow samples (Table 2)   lines in Figure 5. As a result, they may stand for a general formulation for the viscous isotropic behavior of snow according to its porosity through the four parameters (a, b, p, q) given in Table 3. Let us remark that the parameters a and b are close to the ones obtained for metal powders by Geindreau et al. (1999b). However, the exponents corresponding to the snow case are approximately twice bigger, which is linked to a more pronounced dependence on the porosity for very porous materials.
Another interesting feature which can be highlighted in Figure 4 is the fact that the isodissipation curves are closed for all the curves. Their shapes provide information about the ability of an isotropic macroscopic loading to locally activate the ice viscous behavior. Based on the Abouaf formulation (24), the ratio between the maximum isotropic stress S max 1 /3 and the maximum deviatoric stressS 2 max can be expressed for each sample as As snow is always submitted to a mechanical loading which can be decomposed into a deviatoric part and an isotropic part, this 5 ratio provides a measure of the relative contribution of the isotopic part of the mechanical loading in the activation of the ice viscosity. The bigger this ratio, the smaller the activation degree. The evolution of this ratio is plotted in Figure 6 as a function of the snow compacity. The diamond points are computed using the values for f and c presented in Table 2 and the solid line is computed using the two phenomenological fits proposed in (28). The increase in this ratio with snow density highlights the fact that deviatoric fluctuations get smaller under isotropic loading conditions as snow gets denser. In other words, the ice viscosity 10 is more difficult to activate for dense snow than for fresh snow. The divergence of the solid line around 1 corresponds to the limit case of ice where S max 1 becomes infinite as predicted by (3).

Normality rule
The proposed macroscopic modeling is formulated within the framework of associated viscoplasticity. In other words, the flow directionĖ is by construction supposed to be orthogonal to the isodissipation curves (see equation (13)). In the space composed of the two first stress and strain invariants' planes, for an isodissipation curve corresponding to Σ eq , this normality is written In the case of the Abouaf's model (see equation (24)), the theoretical values of the strain flow vectors associated to the seven points in Figure 4 is written The observed difference between theoretical and numerical flow vectors actually results from the slight misfit between the Abouaf models and the numerical points, which is amplified by the radial projection procedure used in order to compute the theoretical flow vectors. Moreover, the validity of normality rule tends to get less accurate as the porosity of the material increases. A similar trend has been already observed in the case of power law fluid flow through porous media (Orgéas et al., 15 2007). Overall, the Abouaf's model presented in section 3.3 provides a satisfactory modeling of the snow viscous behavior on the whole range of investigated densities.

Application to classical laboratory tests
In the case of isotropic snow microstructures, the homogenized constitutive viscous behavior developed in this paper can be summarized as follow: and  where n = 4.5 and A = 1.5 × 10 −3 MPa −4.5 .s −1 account for the ice viscosity and a = 1.5, p = 2.45, b = 8.97 and q = 2.28 account for the snow porosity.
In the following, the mechanical responses of the proposed model are analyzed and compared in the case of classical laboratory tests (Figure 7). In this figure the situation (a) corresponds to the oedometric compression test in which the radial deformation E rr of the snow sample is prevented. The snow mechanical response is then characterized by the relationship between the axial 5 stress Σ zz and the axial strain rateĖ zz . In Figure 7, the situation (b) corresponds to the general triaxial compression test in which the radial stress Σ rr is prescribed and kept constant. From this general setting, two particular cases can be studied: the uniaxial compression test with Σ rr = 0 and the isotropic compression test with Σ zz = Σ rr . In all this section, the classical soil mechanics convention is adopted, i.e. compression stresses are positive, and the elasticity of snow is neglected.

10
As the snow samples are often extracted from the snowpack thanks to hollow cylinders, the oedometric compression test is one of the most convenient mechanical laboratory test to perform on snow. Under the lateral constraint E rr = 0, we have: As a result, the two first strain rate and stress invariants are written 15 In this particular case, from (32) and (33) it can be shown that the lateral constraint E rr = 0 implies that: and consequently, 18 The Cryosphere Discuss., doi:10.5194/tc-2016Discuss., doi:10.5194/tc- -272, 2016 Manuscript under review for journal The Cryosphere anḋ In oedometric experimental tests, the lateral pressure Σ rr is not easily accessible and it is often tempting to neglect this pressure and interpret the oedometric compression test as a uniaxial compression test. The relation (37) can be used to assess the relative importance of the confining pressure with respect to the vertical stress. The evolution of this ratio with respect to 5 the snow density is shown in Figure 8. It should be noticed that this ratio does not depend on the axial strain rateĖ zz . As expected, this ratio increases with increasing the snow density, and tends towards one for φ = 0, due to the incompressiblity of the ice skeleton. The evolution of this ratio is similar to the one measured by Geindreau et al. (1999a) on metallic powders.
In the whole range of snow compacities under consideration (materialized by the grey zone in Figure 8) the lateral pressure represents 30 to 50 % of the vertical stress and cannot be neglected in practice.

Triaxial compression test
During a triaxial test, the cylindric snow sample is submitted simultaneaously to an axial stress Σ zz and a lateral confining pressure Σ rr . Consequently, we have: As a result, the two first strain rate and stress invariants are written In this particular case, from (32) and (33) it can be shown that: and

10
In the case of an uniaxial compression test, Σ rr must be set to 0 in the above equations.
In order to compare the mechanical response of snow under various loadings (uniaxial, oedometric, isotropic and triaxial tests), Figure 9 presents the evolution of the snow densification rate given by E 1 =ρ snow /ρ snow with respect to the snow compacity when constant stresses are applied on the sample. As expected, this figure shows that: whatever the loading, the densification rate strongly decreases with increasing the snow density. In the investigated range, 15 i.e. ρ snow /ρ ice ∈ [0.1, 0.6], the densification rate decreases by 9 orders of magnitude from 10 −1 s −1 to 10 −10 s −1 .
for a given snow density, the loading conditions influence strongly the densification rate. Typically, when Σ zz = 10 kPa the densification rate decreases by nearly one order of magnitude if the confining pressure Σ rr is reduced from 10 kPa (isotropic compression) to 0 kPa (uniaxial compression). On the contrary, when Σ rr = 10 kPa, the densification rate increases by nearly one order of magnitude if the axial stress Σ zz increases from 10 kPa (isotropic compression) to 20 20 kPa (triaxial compression). As expected, this last result shows the increase in the densification rate with the increase in the deviatoric stress (i.e.S 2 ). Even if the lateral confining pressure cannot be neglected during oedometric test as highlighted in Figure 8, the oedometric compression test results in a similar densification rate as the uniaxial compression test for the same axial stress Σ zz . Indeed, the vertical strain rate is lower for an oedometric compression than for a uniaxial one but the geometrical constraint imposed in the oedometric compression test prevents the snow sample from 25 dilating, which is not the case for the uniaxial compression test. Overall the two effects cancel out. Finally, above the classical snow compacity range (ρ snow /ρ ice ≥ 0.6), the densification rate dramatically decreases for the oedometric and isotropic compression tests due to the ice incompressibility. As already underlined in Figure 8 in this limit case, the oedometric compression test is equivalent to the isotropic compression one.
In practice, the strain rate is often imposed on the sample. The Figure 10 presents the evolution of the stress Σ zz versus the snow compacity ρ snow /ρ ice = 1 − φ for two different values of strain ratesĖ zz ∈ {10 −7 ; 10 −5 } s −1 and the different loading 5 conditions (uniaxial, oedometric, isotropic and triaxial tests). This figure suggests the following comments: as expected, for a given strain rate, the stress Σ zz increases with increasing snow density.
for a given strain rate and a given density, the stress Σ zz increases with increasing the lateral pressure Σ rr around the sample.
for a given snow density, the stress Σ zz strongly increases with increasing strain rate, which is in accordance with the the ice viscoplastic behavior is recovered when ρ snow /ρ ice tends towards 1. For a given strain rate, the axial stress for a uniaxial or triaxial compression test tends towards a maximum value. By contrast, due to the incompressibility of ice, the axial stress Σ zz tends towards +∞.
In order to quantitatively compare the predictions of our model against the experimental results of Bartelt and von Moos (2000), we consider a snow of density ρ snow = 255 kg.m −3 (corresponding to ρ snow /ρ ice = 0.27) subjected to a confining pressure 5 of Σ rr = 2.5 kPa and a strain rate of 2.2 × 10 −5 s −1 . In this case, the axial stress predicted by our model is Σ zz = 22.8 kPa, which is consistent with the experimental values obtained by Bartelt and von Moos (2000) around 30 kPa.
Further comparison with Bartelt and von Moos (2000) can be achieved in the case of the uniaxial compression test (Σ rr = 0).
In this case, the axial stress simply reads For a given strain rate, the mechanical response of snow can be compared to the one of ice as in Bartelt and von Moos (2000) by using the following parameter: This parameter compares the axial stress that a given snow sample can transmit (Σ snow zz ) to a rough estimate of this stress given as a fraction of the axial stress transmitted in the case of ice ((1 − φ)Σ ice zz ). In Figure 11, the above theoretical expression of α η (φ) is compared with the experimental fit α η = 0.0028 exp(0.008ρ ice * (1 − φ)) proposed by Bartelt and von Moos (2000). As expected, α η (φ) increases with increasing snow density. By definition, α η (φ) should vary between 0 and 1. We 10 can observe that the theoretical expression of α η (φ) is strictly greater than 1 for ρ snow /ρ ice ∈ [0.8, 1], which is not physically reasonable. This feature results from the independent choices of the parameters a, b, p and q in the fitting procedure used in subsection 4.2. An implicit relation between these parameters could help in order to ensure that α η (φ) remains lower than 1 in the whole compacity range. Nevertheless, in the range of snow densities under consideration (ρ snow /ρ ice ∈ [0.1, 0.6]), α η increases monotonously between roughly 0.1 and 0.6. This prediction is higher than the experimental fit proposed by Bartelt 15 and von Moos (2000) (see Figure 11). But during the experiment, the macroscopic mechanical response probably results from both the viscous deformation of the ice skeleton and the rupture of ice bridges between snow grains after few percents of deformation. This feature is not taken into account in the present model and may explain the observed differences.

Conclusions
Despite the non-linearity of the ice viscous constitutive equation, the homogenization approach introduced by Wautier et al. 20 (2015) is successfully adapted to the numerical homogenization of the snow viscous behavior. By contrast to the elastic case, the macroscopic stress response is not a linear function of the imposed macroscopic strain anymore. As a result, the macroscopic response of snow is investigated in terms of isodissipation curves in the planes of the two first strain rate and stress invariants.
Thanks to a few selected loading paths, an Abouaf's model was fitted onto the numerical results. This formulation seems to be relevant to describe the snow viscoplastic behavior in the whole range of snow density under consideration, provided that the 25 snow microstructure is isotropic. The influence of the snow density, i.e. the porosity, on the viscoplastic response of snow is described through two functions (f (φ) and c(φ), see (28)).
The robustness of this Abouaf formulation is tested for several isotropic snow samples covering the whole range of accessible densities. The fitted models proved to be able to account for the stress and strain rate levels as well as the viscous flow directions.
In particular, the ability of snow to exhibit a viscous behavior even under isotropic strain loading is recovered. The scope of This study Bartelt and von Moos (2000) Figure 11. Comparison between our theoretical prediction for αη (45) and the experimental fit proposed by Bartelt and von Moos (2000).
The compacity range for snow samples under consideration is materialized by the grey zone.
application of the presented unified formulation is quite promising and could help improve the modeling of the densification of the snowpack in avalanche forecasting models.
The proposed homogenization model can be easily used to predict the viscous behavior of snow in classical laboratory tests as illustrated in the last section of this paper. However, the uncertainties made on our model parameters should be quantified through a sensitivity analysis, in order to reckon the ability of our homogenized law for snow viscosity to quantitatively 5 recover the experimental results of Desrues et al. (1980); Bartelt and von Moos (2000); Moos et al. (2003); Scopozza and Bartelt (2003b).
Even though the porosity is known to have a very strong influence on the resulting homogenized properties of snow, it is also acknowledged that the very strong anisotropy of some snow microstructures cannot be neglected. The importance of this anisotropy was quantified within the framework of elasticity (Srivastava et al., 2010(Srivastava et al., , 2016Wautier et al., 2015) but additional 10 work should be carried out in order to extend our homogenized visco-plastic formulation to anisotropic snow types.