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

Antoine Wautier, Christian Geindreau, Frédéric Flin

► To cite this version:

Antoine Wautier, Christian Geindreau, Frédéric Flin. Numerical homogenization of the viscoplastic behavior of snow based on X-ray tomography images. The Cryosphere, 2017, 11, pp.1465-1485. 10.5194/tc-11-1465-2017. insu-03668341

HAL Id: insu-03668341

https://insu.hal.science/insu-03668341

Submitted on 16 May 2022

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

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

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

Antoine Wautier<fnref n="1" />,<fnref n="2" />,a,b,c, Christian Geindreau<fnref n="1" />, and Frédéric Flin<fnref n="2" />

*Correspondence to:* Antoine Wautier (antoine.wautier@irstea.fr)

Received: 24 November 2016 – Discussion started: 5 December 2016
Revised: 16 March 2017 – Accepted: 16 May 2017 – Published: 28 June 2017

**Abstract.** 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 post-processing 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.

## 1 Introduction

Accurately predicting the macroscopic behavior of snow in a wide range of loads, strain rates and temperatures is of a particular interest with respect to avalanche risk forecasting or to structural design of avalanche defense structures. The macroscopic mechanical behavior of snow is, however, strongly influenced by many microscale parameters such as its density

(Mellor, 1974), its microstructure topology (Shapiro et al., 1997) or the mechanical behavior of its ice matrix (elastic, viscoplastic, brittle failure). Disregarding the influence of temperature (Schweizer and Camponovo, 2002), snow exhibits two very different types of mechanical behavior depending on the strain rate (Schulson and Duval, 2009). At high strain rate, large deformations are mainly controlled by bond failures and grain rearrangements, whereas at very low strain rates (typically smaller than $10^{-5} s^{-1}$) snow exhibits a viscoplastic behavior (Narita, 1984; Salm, 1982), which plays an important role in the long-term densification of the snowpack through the microscale deformation of the ice skeleton.

During the last decades, many experiments have been performed to explore the macroscopic behavior of different snow types under various loading conditions and temperatures (Mellor, 1974; Salm, 1982; Desrues et al., 1980; Shapiro et al., 1997; Bartelt and von Moos, 2000; Moos et al., 2003; Scopozza and Bartelt, 2003a). Within the framework of continuum mechanics, several phenomenological models have then been proposed in order to account for these experimental data (Desrues et al., 1980; Scopozza and Bartelt, 2003b; Cresseri and Jommi, 2005; Navarre et al., 2007; Cresseri et al., 2009). The fitted material parameters arising in these models often only 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;

<footnote n="1">¹Université Grenoble Alpes, CNRS, Grenoble INP (Institute of Engineering), 3SR, 38000 Grenoble, France</footnote>

<footnote n="2">²Météo-France – CNRS, CNRM UMR 3589, CEN, 38400 Saint Martin d'Hères, France</footnote>

<footnote>a now at: AgroParisTech-ENGREF, 19 avenue du Maine, 75732 Paris, France</footnote>

<footnote>b now at: Irstea UR RECOVER, 3275 Rte Cézanne, CS 40061, 13182 Aix-en-Provence CEDEX 5, France</footnote>

<footnote>c now at: Université Grenoble Alpes, Irstea, UR ETGR, 2 rue de la Papeterie-BP 76, 38402 St-Martin-d'Hères, France</footnote>
Srivastava et al., 2010; Pinzer et al., 2012; Wang and Baker, 2013; Adams and Walters, 2014; Calonne et al., 2015) good databases of 3-D images for the main snow types reported 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 upscaled 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 numerical studies, discrete element (DEM) and finite element (FEM) methods are the main techniques used to bridge the gap between the topology of the ice skeleton and snow macroscopic mechanical behavior.

- In DEM, the snow skeleton is viewed as an assembly of ice grains interacting with each other at contact points. This method is well suited to model complex interfacial interactions between snow grains, possibly leading to grain rearrangements such as elasto-viscoplastic contact deformation, grain sintering, bond breakage or grain sliding. It has been already used on 3-D idealized assemblies of ice grains to identify microstructural deformation mechanisms and to simulate creep densification processes (Johnson and Hopkins, 2005). Application of DEM directly on X-ray tomography images is, however, not straightforward, since the ice skeleton must be segmented into ice grains approximated with simple geometrical shapes. Recently, DEM simulations taking into account cohesion and friction at the contact between grains have been performed on 3-D assemblies of grains deduced from X-ray tomography and approximated by clumps of spheres (Hagenmuller et al., 2015; Mede et al., 2017). In addition, at the scale of ice grains, the crystalline orientation is known to play a role on the viscous deformation mechanisms at the contact between two grains (secondary creep), as it has been recently reported by Burr et al. (2015b, a, 2017). The full granular structure of snow (grain shape and crystal orientations) can be determined on real 3-D images of snow by X-ray diffraction contrast tomography (DCT) (Rolland du Roscoat et al., 2011) but the application of this technique is not straightforward and very few images are available.

- In FEM, the complex 3-D snow skeleton observed by X-ray microtomography can be meshed without losing any information on the microstructure and different types of mechanical behavior can be assumed for polycrystalline ice. While grain rearrangements cannot be considered without the use of mesh adaptation techniques, the deformation of the ice skeleton is explicitly taken into account, which is of primary importance as long as the topology of the ice skeleton is preserved. In the last decade, many studies addressed in details the case of elasticity (Schneebeli, 2004; Pieritz et al., 2004; Srivas-

tava 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), whereas more complex types of constitutive behavior were hardly ever considered. To the best of our knowledge, no general viscoplastic constitutive equation for snow can be found in the literature. Such an equation would, however, be of particular interest to address the long-term densification of the snowpack under its own weight. During this process, snow deforms at a low strain rate and an important part of the densification results from the viscous deformation (secondary creep) of the ice skeleton. Among the scarce microscale FEM modeling of snow plasticity, Theile et al. (2011) proposed a beam network model based on 3-D images to simulate uniaxial creep of snow, and Chandel et al. (2014) used an elasto-plastic constitutive law for ice in order to determine failure envelopes.

In the following, we formulate a 3-D macroscopic viscoplastic constitutive law for snow by performing FEM simulations on 3-D images. As performed in Wautier et al. (2015) within the framework of elasticity, typical kinematically uniform boundary conditions (KUBC) are applied to 3-D images of snow, and finite element simulations are run in order to link the macroscopic stress response of different snow samples to the imposed strain rates in an incremental form. For each numerical simulation, only small strains are considered to avoid any important microstructure modification. Then, the macroscopic law is generalized to finite deformation problems thanks to the use of a collection of 3-D snow images exhibiting different microstructures and densities. In this upscaling process, the viscous behavior of ice is described by a power law of exponent *n* (secondary creep) as in Theile et al. (2011). Due to the nonlinear behavior under consideration, the homogenization does not provide the complete structure of the macroscopic constitutive equation (Auriault et al., 1992, 2002; Geindreau and Auriault, 1999; Orgéas et al., 2007). Nevertheless, it can be shown that the exponent *n* is preserved at the macroscopic scale and that the macroscopic dissipation power is the volume average of the local one (Sujet, 1993). Using these properties, the macroscopic constitutive equation of snow is formulated within the framework defined by the theory of representation of anisotropic tensor functions (Smith, 1971; Liu, 1982) and by using macroscopic iso-dissipation surfaces (Green, 1972; Abouaf, 1985; Duva and Crow, 1992; Sofronis and McMeeking, 1992; Geindreau et al., 1999b; Storakers et al., 1999; Sanchez et al., 2002; Orgéas et al., 2007).

The paper is organized as follows. In Sect. 2, the numerical homogenization procedure used in Wautier et al. (2015) is recalled and adapted to the study of nonlinear constitutive equations. In Sect. 3, the post-processing procedure used in order to characterize the macroscopic viscous behavior of snow in terms of macroscopic iso-dissipation curves is pre-
sented. These curves might be seen as the equivalent of yield surfaces in plasticity as they characterize the set of stress or strain rate states leading to the same level of mechanical dissipation. Their shapes result from the strong coupling between the microstructure of the ice skeleton and the viscous ice behavior at the microscale. They characterize the 3-D viscoplastic behavior of snow. In Sect. 4, we show that Abouaf models (Abouaf, 1985) are well suited to describe the macroscopic viscous behavior of snow deduced from our numerical simulations. At the end of this section, we propose a macroscopic formulation of the viscoplastic behavior of snow. Finally, in Sect. 5, the mechanical responses of snow for classical experimental tests (uniaxial, oedometric and triaxial compression tests) are modeled thanks to our upscaled law. This illustrates the potential applications of our 3-D homogenized constitutive behavior.

## 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, 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 millimeters 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). In the following, in order to distinguish the two scales of interest, lowercase letters are used for microscopic quantities while uppercase ones are used for their macroscopic counterparts.

Irrespective of the size of the sample considered, the boundary conditions used in a homogenization procedure introduce 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, KUBC with a macroscopic homogeneous strain imposed on the boundary, and periodic boundary conditions (PBCs) with a periodicity condition imposed on the displacement field and the normal stress across the sample boundaries. Although PBCs 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 to assume that the pores are filled by a soft material. In order to avoid

the introduction of such artifacts, KUBC were retained. 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 the 3-D microtomographic images (Step 1) and (ii) defining the kinematic relation $u = E \cdot x$ between the homogeneous macroscopic strain E and the displacement field u on the boundary (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 plugin in Homtools (Lejeunes and Bourgeois, 2011). More details can be found in Wautier et al. (2015). The next two steps (steps 3 and 4) are modified in order to take into account the change in the constitutive modeling of ice.

### 2.1 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 ($\dot{\epsilon}$) is decomposed as the sum of an elastic part ($\dot{\epsilon}_e$) and a viscous part ($\dot{\epsilon}_v$) as

$$ \dot{\epsilon} = \dot{\epsilon}_e + \dot{\epsilon}_v. \qquad (1) $$

The elastic part can be expressed as $\epsilon_e = (C^{\text{ice}})^{-1} : \sigma$, where $C^{\text{ice}}$ is the elastic stiffness tensor, $\sigma$ is the Cauchy stress tensor and “:” the double contraction product. Due to isotropy, $C^{\text{ice}}$ is fully defined by a Young’s modulus E and a Poisson ratio $\nu$. Concerning the viscous part, at low strain rates, the nonlinear mechanical behavior of ice is usually described by a power law (Mellor, 1974; Schulson and Duval, 2009), i.e., the Norton Hoff in 3-D (Lemaitre and Chaboche, 1985). Thus, we have

$$ \dot{\epsilon}_v = \frac{3}{2} A \sigma_{\text{eq}}^{(n-1)} \bar{\sigma}, \qquad (2) $$

where A and n are two material parameters (which usually depend on the temperature), $\bar{\sigma}$ is the deviatoric stress tensor and $\sigma_{\text{eq}}(\sigma)$ is the equivalent stress defined as

$$ \sigma_{\text{eq}}(\sigma) = \sqrt{\frac{3}{2} \bar{\sigma} : \bar{\sigma}} = \sqrt{\frac{3}{2} \bar{s}_2}, \quad \bar{\sigma} = \sigma - \frac{1}{3} \text{Tr}(\sigma) \mathbf{I}, \qquad (3) $$

where $\mathbf{I}$ is the second-order identity tensor, Tr is the trace operator and $\bar{s}_2$ is the second invariant of $\bar{\sigma}$. It can be shown (Lemaitre and Chaboche, 1985) that the viscous strain rate tensor $\dot{\epsilon}_v$ derives from a viscous potential $\omega(\sigma)$ as

$$ \dot{\epsilon}_v = \frac{\partial \omega}{\partial \sigma} = \frac{d\omega}{d\sigma_{\text{eq}}} \frac{\partial \sigma_{\text{eq}}}{\partial \sigma} = \frac{3}{2} \frac{d\omega}{d\sigma_{\text{eq}}} \frac{\bar{\sigma}}{\sigma_{\text{eq}}}, \\ \text{with } \omega(\sigma) = \frac{A}{n+1} \sigma_{\text{eq}}^{(n+1)}. \qquad (4) $$

From Eqs. (2) and (4), one can define the equivalent strain rate $\dot{\epsilon}_{\text{eq}}(\dot{\epsilon}_v)$ as the dual variable of the equivalent stress $\sigma_{\text{eq}}$
<figure><img src="image_4.png" /><figcaption>**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.</figcaption></figure>

such that

$$ \dot{\epsilon}_{eq}(\dot{\epsilon}_v) = \sqrt{\frac{2}{3} \frac{\dot{\epsilon}}{\bar{\epsilon}_v} : \frac{\dot{\epsilon}}{\bar{\epsilon}_v}} = \sqrt{\frac{2}{3} \bar{\epsilon}_2} \quad \text{with} \quad \dot{\epsilon}_{eq} = A\sigma_{eq}^n, \quad (5) $$

where $\bar{\epsilon}_2$ is the second invariant of the deviatoric part of $\dot{\epsilon}_v$. If $p_v$ stands for the volumetric mechanical dissipation, the equivalent stress and the equivalent strain rate verify

$$ p_v = \sigma : \dot{\epsilon}_v = \bar{s}_2 \bar{\epsilon}_2 = \dot{\epsilon}_{eq} \sigma_{eq} = A \sigma_{eq}^{(n+1)}. \quad (6) $$

At the microscopic scale, the viscoplastic deformation of ice is incompressible. Consequently the equivalent stress (or the equivalent strain rate) depends only on the second invariant $\bar{s}_2$ of $\bar{\sigma}$ (or the second invariant $\bar{\epsilon}_2$ of $\dot{\epsilon}_v$).

Overall, the ice matrix is thus modeled as an isotropic elasto-viscoplastic material in the finite element commercial software Abaqus. The values of the material constants used in this constitutive modeling, namely $A, n, E$ and $\nu$, are given in Table 1. It should be underlined that even if the most common values used for $E$ and $n$ are $E = 9$ GPa and $n = 3$, the Young modulus values of ice found in the literature range from 0.2 to 9.5 GPa (Chandel et al., 2014) and the values for $n$ vary between 1.8 and 4.6 under usual loading and temperature conditions (Scapozza and Bartelt, 2003; Schulson and Duval, 2009; Schleef et al., 2014). As a result, different values for $E$ and $n$ are considered in this study to give more insight on their influence on the homogenized viscoplastic behavior of snow.

<table><caption>**Table 1.** Mechanical parameters used in the elasto-viscoplastic modeling of ice implemented in Abaqus.</caption><thead><tr><td>Parameters</td><td>Value</td></tr></thead><tbody><tr><td>A</td><td>1.5 × 10<sup>-3</sup> MPa<sup>-n</sup> s<sup>-1</sup></td></tr><tr><td>n</td><td>2, 3, 4.5</td></tr><tr><td>E</td><td>325 MPa, 9 GPa</td></tr><tr><td>ν</td><td>0.3</td></tr></tbody></table>

## 2.2 Macro-strain paths definition (Step 4)

Given a time-dependent macroscopic strain loading $E(t)$, the KUBC homogenization problem to be solved reads as

$$ \left\{ \begin{array}{ll} \text{div}\sigma = 0 & \text{for } x \in V \\ \mathbf{u} = \mathbf{E}(t) \cdot \mathbf{x} & \text{for } x \in \partial V \\ \boldsymbol{\epsilon} = \frac{1}{2}(\nabla \mathbf{u} + t \nabla \mathbf{u}) & \text{for } x \in V \\ \boldsymbol{\epsilon} = \boldsymbol{\epsilon}_e + \boldsymbol{\epsilon}_v & \text{for } x \in V, \quad (7) \\ \dot{\boldsymbol{\epsilon}}_v = \frac{3}{2} A(x) \sigma_{eq}^{n-1} \bar{\boldsymbol{\sigma}} & \text{for } x \in V \\ \boldsymbol{\sigma} = \frac{E(x)}{1+\nu} \left( \boldsymbol{\epsilon}_e + \frac{\nu}{1-2\nu} \text{Tr}(\boldsymbol{\epsilon}_e) \mathbf{1} \right) & \text{for } x \in V \end{array} \right. $$

where $V$ stands for the domain occupied by the whole snow sample and $\partial V$ its boundary. The spatial heterogeneity of the
mechanical properties of snow is captured thanks to the functions $A(\mathbf{x})$ and $E(\mathbf{x})$ defined as

$$ 
\begin{aligned}
A(\mathbf{x}) &= \begin{cases} A & \text{if } \mathbf{x} \in V_i \\ 0 & \text{if } \mathbf{x} \in V \setminus V_i \end{cases} \quad \text{and} \\
E(\mathbf{x}) &= \begin{cases} E & \text{if } \mathbf{x} \in V_i \\ 0 & \text{if } \mathbf{x} \in V \setminus V_i \end{cases},
\end{aligned} 
\quad (8) 
$$

where $V_i \subset V$ is the domain occupied by the ice matrix. Similarly to the elastic case (Wautier et al., 2015), the macroscopic stress tensor $\Sigma$ is deduced from the knowledge of its microscopic counterpart thanks to the volume averaging:

$$ \Sigma = \frac{1}{|V|} \int_V \sigma dV = \langle \sigma \rangle. \quad (9) $$

As a result, for a given macroscopic strain loading $E(t)$, the macroscopic stress response $\Sigma(t)$ is recovered. The implicit function linking these two second-order tensors characterizes the homogeneous behavior of the snow sample considered and can be put in the form:

$$ \dot{\mathbf{E}} = \dot{\mathbf{E}}_e + \dot{\mathbf{E}}_v = \mathcal{F}(\mathbf{\Sigma}), \quad (10) $$

where $\dot{\mathbf{E}}_e$ is the macroscopic elastic strain rate tensor and $\dot{\mathbf{E}}_v$ is the macroscopic viscous strain rate tensor. The elastic part can be expressed as $\mathbf{E}_e = (\mathbf{C}^{\text{hom}})^{-1} : \mathbf{\Sigma}$, where $\mathbf{C}^{\text{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 3-D images. A single simulation is required if the snow microstructure is isotropic. By contrast, the homogenization of the viscoplastic 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; Orgéas et al., 2007). Indeed, the following can be shown:

- The homogeneity of degree $n$ of the microscopic viscous constitutive Eq. (2) is preserved in the homogenization process. In other words, the macroscopic viscous strain rate $\dot{\mathbf{E}}_v$ is a homogeneous function of degree $n$ of the macroscopic stress $\mathbf{\Sigma}$, and the macroscopic volumetric mechanical dissipation $\mathcal{P}_v = \dot{\mathbf{E}}_v : \mathbf{\Sigma}$ is an homogeneous function of degree $n+1$ of $\mathbf{\Sigma}$.

$$ \forall \lambda \in \mathbb{R} \quad \left\{ \begin{aligned} \dot{\mathbf{E}}_v(\lambda \mathbf{\Sigma}) &= \lambda^n \dot{\mathbf{E}}_v(\mathbf{\Sigma}) \\ \mathcal{P}_v(\lambda \mathbf{\Sigma}) &= \lambda \mathbf{\Sigma} : \dot{\mathbf{E}}_v(\lambda \mathbf{\Sigma}) = \lambda^{n+1} \mathcal{P}_v(\mathbf{\Sigma}). \end{aligned} \right. \quad (11) $$

As a result, the choice in the macroscopic strain rate $\dot{\mathbf{E}}_v$ can be reduced to the unit sphere in the second-order tensor space, i.e., to strain rate tensors of norm $\sqrt{\dot{\mathbf{E}}_v} : \dot{\mathbf{E}}_v = 1$.

- The macroscopic dissipation potential $\Omega(\mathbf{\Sigma})$ is the volume-average of the local dissipation potential $\omega$:

$$ \Omega(\mathbf{\Sigma}) = \frac{1}{|V|} \int_V \omega(\mathbf{\sigma}) dV = \langle \omega(\mathbf{\sigma}) \rangle. \quad (12) $$

Consequently, as at the microscopic scale (see Eq. 4), we have

$$ \begin{gathered} \dot{\mathbf{E}}_v = \frac{\partial \Omega}{\partial \mathbf{\Sigma}} = \frac{d\Omega}{d\Sigma_{\text{eq}}} \frac{\partial \Sigma_{\text{eq}}}{\partial \mathbf{\Sigma}}, \\ \text{with } \Omega(\mathbf{\Sigma}) = \frac{A}{n+1} \Sigma_{\text{eq}}^{n+1}, \end{gathered} \quad (13) $$

where $\Sigma_{\text{eq}}(\mathbf{\Sigma})$ is the macroscopic equivalent stress. $\Sigma_{\text{eq}}(\mathbf{\Sigma})$ verifies

$$ \mathcal{P}_v = \mathbf{\Sigma} : \dot{\mathbf{E}}_v = \dot{\mathbf{E}}_{\text{eq}} \mathbf{\Sigma}_{\text{eq}} = A \mathbf{\Sigma}_{\text{eq}}^{n+1}, \quad (14) $$

with $\dot{\mathbf{E}}_{\text{eq}}(\dot{\mathbf{E}}_v)$ the macroscopic equivalent strain rate defined by duality.

As a result, the macroscopic viscoplastic law (Eq. 13) is perfectly defined if the macroscopic equivalent stress $\Sigma_{\text{eq}}$ is known. Equation (14) shows that this macroscopic equivalent stress $\Sigma_{\text{eq}}$ can be fitted on isomechanical dissipation surfaces in the space associated with $\mathbf{\Sigma}$. Let us remark that the shape and size of such isomechanical dissipation surfaces result from the strong coupling between the microstructure and the nonlinear behavior of the ice under consideration. The relation (Eq. 14) also shows the equivalent stress can be obtained whatever the chosen $A$ value. In the case of general anisotropy, the form of $\Sigma_{\text{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 $\mathbf{\Sigma}$ as a consequence of definition (Eq. 4). This condition may be expressed as $\Omega(\mathbf{\Sigma}) = \Omega(-\mathbf{\Sigma})$. Finally, let us remark that by definition (see Eqs. 13 and 14) the macroscopic strain rate $\dot{\mathbf{E}}_v$ is normal to isomechanical 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, for a given value of $n$, it can be shown (Abouaf, 1985; Geindreau et al., 1999b; Danas et al., 2008) that the macroscopic equivalent stress is written

$$ \Sigma_{\text{eq}}(\mathbf{\Sigma}) = \Sigma_{\text{eq}}(S_1, \bar{S}_2, \bar{S}_3, \phi), \quad (15) $$
where $\phi$ is the snow porosity and $(S_1, \bar{S}_2, \bar{S}_3)$ are the three invariants of the macroscopic stress tensor $\Sigma$ defined as

$$ S_1 = \text{Tr}(\Sigma), \quad \bar{S}_2 = \sqrt{\bar{\Sigma} : \bar{\Sigma}}, \quad \bar{S}_3 = \det(\bar{\Sigma}), \\ \text{with } \bar{\Sigma} = \Sigma - \frac{S_1}{3} \mathbf{I}. \qquad (16) $$

Similarly, the macroscopic equivalent strain rate $\dot{E}_{eq}$ takes the form

$$ \dot{E}_{eq}(\dot{E}_{v}) = \dot{E}_{eq}(E_{1}, \bar{E}_{2}, \bar{E}_{3}, \phi), \qquad (17) $$

where $(E_1, \bar{E}_2, \bar{E}_3)$ are the invariants of the strain rate tensor $\dot{E}_v$ defined as

$$ E_1 = \text{Tr}(\dot{E}_v), \quad \bar{E}_2 = \sqrt{\dot{\bar{E}}_v : \dot{\bar{E}}_v}, \quad \bar{E}_3 = \det(\dot{\bar{E}}_v), \\ \text{with } \dot{\bar{E}}_v = \dot{E}_v - \frac{E_1}{3} \mathbf{I}. \qquad (18) $$

In contrast with the microscopic scale (see Eq. 3), the macroscopic equivalent stress (Eq. 15) depends on the three invariants $(S_1, \bar{S}_2, \bar{S}_3)$ of the macroscopic stress tensor $\Sigma$. 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 $(\bar{S}_2)$. The third invariant $\bar{S}_3$ characterizes the loading type and is linked to the Lode angle $\vartheta$ in the stress space (Lemaitre and Chaboche, 1985; Danas et al., 2008):

$$ \cos(3\vartheta) = \frac{27}{2} \frac{\bar{S}_3}{\Sigma_{eq}^3}. \qquad (19) $$

As a first approximation, it seems reasonable to assume that the influence of the third invariant $\bar{S}_3$ is negligible (Green, 1972; Abouaf, 1985; Geindreau et al., 1999b; Fritzen et al., 2012). Consequently, the macroscopic volumetric mechanical dissipation $\mathcal{P}_v$ depends on the first and second stress and strain invariants and not only on the second ones as at the microscale 6.

$$ \mathcal{P}_v = \Sigma : \dot{\mathbf{E}}_v = E_{eq} \Sigma_{eq} = A \Sigma_{eq}^{n+1} = \frac{1}{3} E_1 S_1 + \bar{E}_2 \bar{S}_2 \qquad (20) $$

The relation (Eq. 20) shows that, for a given snow porosity, the equivalent macroscopic stress $\Sigma_{eq}$ can be fitted on iso-volumetric mechanical dissipation curves in the plane $(S_1/3, \bar{S}_2)$. These iso-dissipation curves can be obtained by plotting the values $(S_1/3, \bar{S}_2)$ corresponding to different loading conditions defined by $(E_1, \bar{E}_2)$. Therefore, the choice was made to run numerical simulations for seven diagonal strain rate tensors defined such that the loading direction in the plane $(E_1, \bar{E}_2)$ varies from 0 to 90°. More explicitly, $\dot{\mathbf{E}}$ applied on the sample is taken as

$$ \dot{\mathbf{E}} = \dot{\mathbf{E}}_{\text{ref}} \begin{pmatrix} 1 & 0 & 0 \\ 0 & \eta & 0 \\ 0 & 0 & \eta \end{pmatrix}, \qquad (21) $$

with $\dot{\mathbf{E}}_{\text{ref}} = 10^{-7} \text{s}^{-1}$ and $\eta \le 1$ such that $\frac{\bar{E}_2}{\bar{E}_1} = \sqrt{\frac{2}{3} \frac{1-\eta}{1+2\eta}} = \tan\theta, \theta \in \{0^\circ, 9^\circ, 18^\circ, 30^\circ, 45^\circ, 65^\circ, 90^\circ\}$.

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 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 $(\ell_1, \ell_2, \ell_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 3-D 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 viscoplasticity

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, \bar{E}_2)$ given by the Eq. (21). The strain rate is applied on each sample during less than 40 000 s, corresponding to a volumetric strain smaller than 1.2%. The overall viscous behavior of the snow samples is deduced thanks to a post-processing procedure consisting in the three steps described in this section (steps a to c) and summarized in Fig. 2.

#### 3.1 Extracting the viscous response (Step a)

Because snow is locally modeled as an elasto-viscoplastic material in Abaqus (see Sect. 2.1), the macroscopic time response $\Sigma(t) = \frac{1}{|V|} \int_V \sigma(x) dV$ deduced from numerical simulations does not involve only viscoplasticity. In the case where the ice viscosity is activated everywhere in the ice skeleton, the macroscopic stress $\Sigma(t)$ corresponding to a constant strain rate should stabilize around a constant value according to Eq. (2). However, due to the complex geometry of the ice skeleton, the ice viscosity is not uniformly activated and the time response of $\Sigma(t)$ is influenced by the ice elastic behavior even in the long term. Based on the material parameters $A, n$ and $E$, and the typical imposed strain rate $\dot{\mathbf{E}}_{\text{ref}}$, a characteristic time $\tau$ can be introduced as the ratio between the ice viscosity $\eta(\dot{\mathbf{E}}_{\text{ref}}) = (\dot{\mathbf{E}}_{\text{ref}}/A)^{1/n}/\dot{\mathbf{E}}_{\text{ref}}$ and the Young modulus $E$:

$$ \tau = \frac{\eta(\dot{\mathbf{E}}_{\text{ref}})}{E} = \frac{1}{E} \left( A^{-\frac{1}{n}} \dot{\mathbf{E}}_{\text{ref}}^{\frac{1-n}{n}} \right). \qquad (22) $$

The typical stress response of a snow sample under a constant given strain rate versus the dimensionless time $(t/\tau)$ is illus-
<table><caption>**Table 2.** Names and principal characteristics of the six snow images of Wautier et al. (2015) used in this study.</caption><thead><tr><th>Sample name</th><th>Snow type</th><th>Dim (px)</th><th>Dim (mm)</th><th>Resolution (µm px<sup>-1</sup>)</th><th>Snow density (kg m<sup>-3</sup>)</th><th>Porosity</th><th>Correlation lengths (ℓ<sub>1</sub>, ℓ<sub>2</sub>, ℓ<sub>3</sub>) (µm)</th></tr></thead><tbody><tr><td>PP_123kg_600</td><td>PP</td><td>600</td><td>2.95</td><td>4.91</td><td>123.31</td><td>0.87</td><td>(64, 64, 65)</td></tr><tr><td>RG_172kg_600</td><td>RG</td><td>600</td><td>2.95</td><td>4.91</td><td>172.74</td><td>0.81</td><td>(92, 94, 97)</td></tr><tr><td>RG_256kg_512</td><td>RG</td><td>512</td><td>2.51</td><td>4.91</td><td>256.28</td><td>0.72</td><td>(113, 111, 110)</td></tr><tr><td>RG_1600</td><td>RG</td><td>600</td><td>4.46</td><td>7.43</td><td>330.13</td><td>0.64</td><td>(117, 111, 108)</td></tr><tr><td>RG_430kg_651</td><td>RG</td><td>651</td><td>5.61</td><td>8.61</td><td>430.59</td><td>0.53</td><td>(83, 82, 81)</td></tr><tr><td>MF_522kg_542</td><td>MF</td><td>542</td><td>5.42</td><td>10.00</td><td>522.31</td><td>0.43</td><td>(138, 134, 133)</td></tr></tbody></table>

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

**Step a:** - Compute the mechanical response ($S_1(t), \bar{S}_2(t)$) for a given constant strain rate $\dot{\bar{E}}$

- Determine the viscous response ($S_1^v, \bar{S}_2^v$) associated with the strain state ($E_1, \bar{E}_2$)

**Step b:** - For each loading given in Eq. (21), compute the associated viscous dissipated power with Eq. (23)

- Rescale the numerical points on an isodissipation curve thanks to Eq. (24)

**Step c:** - Fit an Abouaf model thanks to Eqs. (25) and (28)

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

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

<figure><img src="image_8.png" /><figcaption>**Figure 2.** Three-step post-processing procedure used in order to formulate a homogenized viscous constitutive equation.</figcaption></figure>

trated in Fig. 3 for the snow sample RG_1600 (Table 2) for different values of *n* and *E*. In all the cases, the mechanical response is characterized by a transient regime driven by the elastic properties followed by a permanent regime dominated by the viscoplastic behavior. As illustrated by the comparison between the cases ($n, E) = (4.5, 325$ MPa) and ($n, E) = (4.5, 9$ GPa), the responses $S_1(t/\tau)$ and $\bar{S}_2(t/\tau)$ are independent of the Young modulus value chosen. On the contrary, the mechanical response is influenced by the *n* value.

As a result, for a given value of *n* and whatever the Young modulus, the viscoplastic behavior of snow can be characterized by computing the intersection point between initial and final asymptotes of the curves $S_1(t/\tau)$ and $\bar{S}_2(t/\tau)$ in Fig. 3. The obtained values for the two stress invariants are noted $S_1^v$ and $\bar{S}_2^v$ and are systematically used in the rest of this paper as

the snow viscous homogeneous response to a given imposed constant macroscopic strain rate. The top graph in Fig. 2 illustrates this step of the post-processing procedure.

## 3.2 Computing iso-dissipation curves (Step b)

For a given snow sample of porosity $\phi$ and for each applied loading path (Eq. 21), the macroscopic volumetric mechanical dissipation $P_v$ (Eq. 20) is computed as

$$P_v = \frac{1}{3} E_1 S_1^v + \bar{E}_2 \bar{S}_2^v, \qquad (23)$$

where $S_1^v$ and $\bar{S}_2^v$ are the characteristic stress invariants obtained in the first step of the post-processing process (Fig. 2). Each loading path leads to different values of $P_v$. However,
<figure><img src="image_9.png" /><figcaption>Figure 3. Imposed strain rate (a) and stress response (b) versus dimensionless time (t/τ) for the sample RG_1600 (see Table 2). The loading strain rate is characterized by θ = 65° in Eq. (21). Two Young's moduli and three values of n are considered.</figcaption></figure>

iso-mechanical dissipation points in the plane (S₁/3, S̄₂) can be recovered thanks to the homogeneity property (Eq. 11). Given an arbitrary value of P̄v = 1 Pa s⁻¹, the corresponding macroscopic strain and stress invariants are computed as

$$ (E_1^\circ, \bar{E}_2^\circ) = \left( \frac{P_v^\circ}{P_v} \right)^{n/(n+1)} (E_1, \bar{E}_2), \quad \text{and} \\ (S_1^\circ, \bar{S}_2^\circ) = \left( \frac{P_v^\circ}{P_v} \right)^{1/(n+1)} (S_1^v, \bar{S}_2^v). \qquad (24) $$

Thanks to this rescaling, the seven homogenization tests (Eq. 21) enable the description of an iso-dissipation curve in the plane of the stress invariants (S₁/3, S̄₂) as illustrated in Fig. 2 (Step b). For each point (S̄₁/3, S̄₂) on this graph, the associated flow vector (Ē₁, Ē₂) is plotted. The viscous dissipated power is thus simply equivalent to the scalar product

$$ (S_1^\circ/3, \bar{S}_2^\circ) \cdot (E_1^\circ, \bar{E}_2^\circ). $$

### 3.3 Abouaf's model (Step c)

Within the framework presented in Sect. 2.2, Abouaf (1985) suggested using the macroscopic equivalent stress initially proposed by Green (1972) to describe the viscoplastic behavior of metal powders at high temperatures. This macroscopic

equivalent stress Σeq(Σ), is written

$$ \Sigma_{\text{eq}}(\boldsymbol{\Sigma}) = \Sigma_{\text{eq}}(S_1, \bar{S}_2, \phi) = \sqrt{f(\phi) S_1^2 + \frac{3}{2} c(\phi) \bar{S}_2^2}, \quad (25) $$

where f(φ) and c(φ) are two material functions which depend on snow porosity only for a given exponent n of the constitutive Eq. (2). When φ = 0, we have f(φ) = 0 and c(φ) = 1 in order to recover the equivalent viscous stress of the ice matrix (Eq. 3): Σeq(Σ, φ = 0) = σeq(σ). From the definition of the viscous strain in Eq. (13) together with the previous definition of the equivalent stress in Eq. (25), it can be shown that

$$ \dot{\mathbf{E}}_v = A \Sigma_{\text{eq}}^{n-1} \left( f(\phi) S_1 \mathbf{I} + \frac{3}{2} c(\phi) \bar{\mathbf{\Sigma}} \right). \qquad (26) $$

As a result, the corresponding macroscopic equivalent strain rate introduced in Eq. (20) reads

$$ \dot{\mathbf{E}}_{\text{eq}}(\dot{\mathbf{E}}_v) = \dot{\mathbf{E}}_{\text{eq}}(E_1, \bar{E}_2, \phi) = \sqrt{\frac{E_1^2}{9f(\phi)} + \frac{2}{3} \frac{\bar{E}_2^2}{c(\phi)}}. \qquad (27) $$

For a given porosity φ, the combination of Eqs. (25) and (20) provides an implicit definition of f(φ) and c(φ) such that for
all ($S_1, \bar{S}_2$)

$$ \Sigma_{\text{eq}}(S_1, \bar{S}_2, \phi) = \sqrt{f(\phi) S_1^2 + \frac{3}{2} c(\phi) \bar{S}_2^2} \\ = \left(\frac{P_v^\circ}{A}\right)^{1/(n+1)}. \qquad (28) $$

In the present work, optimal values for $f(\phi)$ and $c(\phi)$ were obtained in the range $\phi \in [0.43, 0.87]$ by minimizing the quadratic error between the model (Eq. 28) and the numerical points ($S_1^\circ(\theta)/3, \bar{S}_2^\circ(\theta)$).

## 4 Results and discussion

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 Fig. 4, the seven points ($S_1^\circ/3, \bar{S}_2^\circ$) corresponding to the strain rates of Eq. (21) are represented for these six snow samples in the plane of the two first stress invariants for $n = 4.5$. Similar results have been obtained for the other values of $n$, as shown in Figure 5. The corresponding strain flow vectors ($E_1^\circ, \bar{E}_2^\circ$) are shown by solid arrows and iso-dissipation curves corresponding to fitted Abouaf models are represented by solid lines. Optimal values for $f$ and $c$ obtained for each snow type are presented in Fig. 4 and reported in Table 3 for $n \in \{2, 3, 4.5\}$. It should be underlined here that each iso-dissipation curve is typical of a given snow characterized by its density and thus each curve is also an iso-density curve.

### 4.1 Iso-dissipation curves for various snow samples

The comparison between simulated points and fitted Abouaf models shows overall good agreement even if a systematic slight deviation is observed for the highest $S_1$ values. It should also be noticed that the stress state corresponding to an isotropic strain rate ($\bar{E}_2 = 0$) is not completely isotropic ($\bar{S}_2 \neq 0$). This feature cannot be captured by Abouaf models, which assumes a perfect isotropic behavior of the material. 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 Sect. 2.2, the viscous behavior of snow should be insensitive to the sign of $\Sigma$ as the ice matrix behaves exactly the same in tension and in compression. In the stress space ($S_1/3, \bar{S}_2$), this results in the symmetry of the iso-dissipation curves with respect to the axis $S_1/3 = 0$. Provided that the iso-dissipation curves are smooth, their tangent for $S_1/3 = 0$ is horizontal, which is fulfilled in Fig. 4. It must be mentioned that when snow is subjected to large strain levels, geometrical effects will introduce nonlinear 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.

The overall viscoplastic response of snow is of course sensitive to the $n$ exponent of Norton Hoff's law (Eq. 2) used for the ice (see Sect. 3.1). As for example in Fig. 5, the influence of $n$ onto the iso-dissipation curves is shown for the snow sample MF_522. Similar results have been obtained on the other samples. As expected, for a given value $P_v^\circ$, the size of the iso-dissipation curves increases with $n$ (since the ice viscosity $\eta(\dot{E}_{\text{ref}})$ increases) but their shape remains unchanged. They can be deduced from each other by simple dilation.

### 4.2 Density dependence of the iso-dissipation curves

As snow density increases, iso-dissipation curves tend to expand, and conversely flow vectors tend to get smaller. In terms of physics, this means that the denser the snow is, the smaller the applied strain rate is 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 snow viscous behavior is fully described by the evolution of $f(\phi)$ and $c(\phi)$ with respect to snow relative density ($\rho_{\text{snow}}/\rho_{\text{ice}} = 1 - \phi$) represented in Fig. 6 for $n \in \{2, 3, 4.5\}$. As snow density increases in Fig. 6, parameter values decrease, which is consistent with the implicit definition of the iso-dissipation curves in Eq. (25). For a given equivalent stress $\Sigma_{\text{eq}}$, higher values for $f$ and $c$ will result in lower stress invariants $S_1$ and $\bar{S}_2$ as observed in Fig. 4. Concerning the influence of $n$, the observed increase in $f$ and $c$ is consistent with the dilation of the iso-dissipation curves observed in Fig. 5.

Different expressions of the material functions $f(\phi)$ and $c(\phi)$ have been proposed in the past based on experimental data on metal powders (Abouaf, 1985; Abouaf and Chenot, 1988; Geindreau et al., 1999b), micromechanical modeling (cell model; Duva and Crow, 1992) or numerical simulations on simple microstructures (Sofronis and McMeeking, 1992). These functions have been identified in a restricted range of porosity (dense materials with $\phi < 0.4$). We propose to fit our numerical results using the expressions proposed by Geindreau et al. (1999b). In order to account for the change in the porosity range between metal powders and snow, the relative density limit value of 0.57 is set equal to zero. As a result, the proposed functions are written

$$ \left\{ \begin{array}{l} f(\phi) = a \left( \frac{\phi}{1-\phi} \right)^p \\ c(\phi) = 1 + b \left( \frac{\phi}{1-\phi} \right)^q, \quad (a, b, p, q) \in \mathbb{R}^4. \end{array} \right. \qquad (29) $$

The above fits respect the theoretical values $f(0) = 0$ and $c(0) = 1$ already mentioned in Sect. 3.3. For highly porous snow ($\phi \to 1$), an infinitely small stress level would be needed in order to produce a high viscous dissipation. This
<figure><img src="image_10.png" /></figure>

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

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

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

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

<figure><img src="image_15.png" /><figcaption>**Figure 4.** Iso-dissipation curves in the plane of the stress invariants ($S_1, \bar{S}_2$) corresponding to $n = 4.5$ and to an arbitrary dissipated power $P_V^°$ for the six snow samples of Table 2 in increasing density order. The associated strain flow vectors ($E_1^°, \bar{E}_2^°$) are represented by solid arrows. Abouaf models are fitted to the numerical points (solid lines) and theoretical values of strain flow vectors are shown (dashed arrows). Parameters $f$ and $c$ of the fits are shown on the graphs and summarized in Table 3.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>**Figure 5.** Influence of the exponent *n* onto the iso-dissipation curves ($P_v = P_v^°$) for the particular snow sample *MF_522*. The associated strain flow vectors ($E_1^°, E_2^°$) are represented by solid arrows. Abouaf models are fitted to the numerical points (solid lines) and theoretical values of strain flow vectors are shown (dashed arrows).</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>**Figure 6.** Evolution of the Abouaf coefficients *f* and *c* (numerical results as diamond points; functions corresponding to Eq. 29 as solid lines) with respect to snow relative density for different *n* values.</figcaption></figure>
<table><caption>**Table 3.** Optimal values for the parameters *f* and *c* of the Abouaf equivalent stress (Eq. 25) for three *n* values.</caption><thead><tr><th rowspan="2">Sample name</th><th rowspan="2">Porosity</th><th colspan="2">n = 2</th><th colspan="2">n = 3</th><th colspan="2">n = 4.5</th></tr><tr><th>f</th><th>c</th><th>f</th><th>c</th><th>f</th><th>c</th></tr></thead><tbody><tr><td>PP_123kg_600</td><td>0.87</td><td>36.0</td><td>150</td><td>79.7</td><td>336</td><td>146</td><td>628</td></tr><tr><td>RG_172kg_600</td><td>0.81</td><td>16.3</td><td>75.7</td><td>33.5</td><td>156</td><td>58.3</td><td>277.4</td></tr><tr><td>RG_256kg_512</td><td>0.72</td><td>4.05</td><td>20.5</td><td>6.98</td><td>34.7</td><td>10.5</td><td>52.3</td></tr><tr><td>RG_1600</td><td>0.64</td><td>2.07</td><td>11.0</td><td>3.32</td><td>17.0</td><td>4.70</td><td>24.0</td></tr><tr><td>RG_430kg_651</td><td>0.53</td><td>0.915</td><td>6.38</td><td>1.40</td><td>9.12</td><td>1.89</td><td>12.1</td></tr><tr><td>MF_522kg_542</td><td>0.43</td><td>0.354</td><td>3.32</td><td>0.503</td><td>4.26</td><td>0.630</td><td>5.07</td></tr></tbody></table>

is consistent with the infinitely high values for *f* and *c* proposed by the above functions (Eq. 29). These functions allow a good description of the numerical points resulting from the homogenization of the six snow samples (Table 3) and are represented by solid lines in Fig. 6 for $n \in \{2, 3, 4.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 *n*-dependent parameters (*a*, *b*, *p*, *q*) given in Table 4. For the sake of illustration, the evolution of the coefficients *a*, *p*, *b* and *q* with respect to *n* is shown in Fig. 7. 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 Fig. 4 is the fact that the iso-dissipation curves are closed for all the snow samples under consideration. This contrasts with the shape of ice iso-dissipation curves in the plane ($S_1/3, \bar{S}_2$) which are represented by horizontal straight lines. Indeed, the corresponding viscous power does not depend on $S_1$ since the constitutive equation (Eq. 2) for ice only involves the deviatoric stress $\bar{\sigma}$. The ability of snow to dissipate significant energy even under an isotropic loading ($\bar{S}_2 = 0$) is due to its porosity. Indeed, even under this type of macroscopic loading conditions, some regions of the microstructure experience a nonzero deviatoric loading, which activates locally the viscous ice behavior. Even if this deviatoric loading vanishes on average, the mean viscous dissipated power only piles up, which results in closed iso-dissipation curves. Their shapes provide information about the ability of an isotropic macroscopic loading to locally activate the viscous ice behavior. Based on the Abouaf formulation (Eq. 25), the ratio between the maximum isotropic stress $S_1^{\max}/3$ and the maximum deviatoric stress $\bar{S}_2^{\max}$ can be expressed for each sample as

$$ \frac{(S_1^{\max}/3)}{\bar{S}_2^{\max}} = \sqrt{\frac{c}{6f}}. \qquad (30) $$

As snow is always submitted to a mechanical loading which can be decomposed into a deviatoric part and an isotropic part, this ratio provides a measure of the relative contribu-

<table><caption>**Table 4.** Optimal parameters chosen for the expressions Eq. (29) for different *n* values.</caption><thead><tr><th></th><th>n = 2</th><th>n = 3</th><th>n = 4.5</th></tr></thead><tbody><tr><th>a</th><td>0.68</td><td>1.0</td><td>1.5</td></tr><tr><th>p</th><td>2.1</td><td>2.3</td><td>2.5</td></tr><tr><th>b</th><td>4.0</td><td>6.1</td><td>8.9</td></tr><tr><th>q</th><td>2.0</td><td>2.2</td><td>2.3</td></tr></tbody></table>

<figure><img src="image_18.png" /><figcaption>**Figure 7.** Evolution of the fitted coefficients *a*, *p*, *b* and *q* with respect to the exponent *n* of the Norton Hoff ice constitutive behavior.</figcaption></figure>

tion of the isotropic part of the mechanical loading in the activation of the ice viscosity. The bigger this ratio is, the smaller the activation degree. The evolution of this ratio is plotted in Fig. 8 as a function of snow relative density for $n \in \{2, 3, 4.5\}$. The diamond points are computed using the values for *f* and *c* presented in Table 3 and the solid line is computed using the two functions (Eq. 29) with the parameters presented in Table 4. 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 is more difficult to activate for dense snow than for fresh snow under isotropic loading conditions. The divergence of the solid line around 1 corresponds to the limit case of ice where $S_1^{\max}$ becomes infi-
<figure><img src="image_19.png" /><figcaption>Figure 8. Aspect ratios ($\bar{S}_1^{\max}/3$)/$S_2^{\max}$ of the Abouaf models with respect to snow relative density. Numerical results are represented with diamond points, results computed using functions (Eq. 29) with solid lines.</figcaption></figure>

nite as predicted by Eq. (3). In Fig. 8, the dependence on the n value is very limited, which highlights the fact that n does not have any influence on the shape of the iso-dissipation curves but only on the stress level leading to a given iso-dissipation as already seen in Fig. 5.

## 4.3 Normality rule

The proposed macroscopic modeling is formulated within the framework of associated viscoplasticity. In other words, the flow direction $\dot{E}$ is by construction supposed to be orthogonal to the iso-dissipation curves (see Eq. 13). In the space composed of the two first stress and strain invariant planes, for an iso-dissipation curve corresponding to $\Sigma_{eq}$, this normality is written

$$ \left( \frac{E_1}{\bar{E}_2} \right) \propto \left( \frac{\frac{3\partial\Sigma_{eq}}{\partial S_1}}{\frac{\partial\Sigma_{eq}}{\partial \bar{S}_2}} \right). \quad (31) $$

In the case of Abouaf models (see Eq. 25), theoretical strain flow vectors associated with the seven points in Fig. 4 are given as

$$ \left( \frac{E_1}{\bar{E}_2} \right) = A \Sigma_{eq}^{n-1} \left( \frac{3f(\phi)S_1}{\frac{3}{2}c(\phi)\bar{S}_2} \right). \quad (32) $$

In Figs. 4 and 5, theoretical strain flow vectors are represented by dashed arrows for the radial projections of the numerical points on Abouaf fits. The overall comparison

with their numerical counterparts represented by solid arrows is quite satisfactory, especially for the densest snows. However, concerning the flow direction, the Abouaf model tends to overpredict the strain deviatoric component for high-deviatoric stresses $\bar{S}_2$ and to underpredict the strain deviatoric component for low-deviatoric stresses except for the direction $\theta = 0^\circ$. In terms of magnitude, the Abouaf model tends to underpredict the intensity of the flow for deviatoric loading.

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., 2007). Overall, the Abouaf model presented in Sect. 3.3 provides a satisfactory modeling of snow viscous behavior on the whole range of investigated densities.

## 5 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

$$ \dot{E}_v = \dot{E}_{eq} \frac{\partial \Sigma_{eq}}{\partial \Sigma} = A \Sigma_{eq}^{n-1} \left( f(\phi) S_1 \mathbf{I} + \frac{3}{2} c(\phi) \bar{\Sigma} \right), \quad (33) $$
with

$$
\begin{align}
\Sigma_{\text{eq}}(S_1, \bar{S}_2, \phi) &= \sqrt{f(\phi) S_1^2 + \frac{3}{2} c(\phi) \bar{S}_2^2}, \\
\dot{E}_{\text{eq}}(E_1, \bar{E}_2, \phi) &= \sqrt{\frac{E_1^2}{9f(\phi)} + \frac{2}{3} \frac{\bar{E}_2^2}{c(\phi)}}, \tag{34}
\end{align}
$$

and

$$
f(\phi) = a \left( \frac{\phi}{1-\phi} \right)^p, \quad c(\phi) = 1 + b \left( \frac{\phi}{1-\phi} \right)^q, \qquad (35)
$$

where *n* and *A* account for the ice viscosity (Table 1) and *a*, *p*, *b* and *q* account for snow porosity (Table 4).

In the following, the mechanical responses of the proposed model are analyzed and compared in the case of classical laboratory tests. In Fig. 9, situation (a) corresponds to an oedometric compression test in which the radial deformation $E_{rr}$ of the snow sample is prevented. Snow mechanical response is then characterized by the relationship between the axial stress $\Sigma_{zz}$ and the axial strain rate $\dot{E}_{zz}$. Situation (b) corresponds to a general triaxial compression test in which the radial stress $\Sigma_{rr}$ is prescribed and kept constant. From this general setting, two particular cases can be studied: a uniaxial compression test when $\Sigma_{rr} = 0$ and an isotropic compression test when $\Sigma_{zz} = \Sigma_{rr}$. In this section, classical soil mechanics convention is adopted, i.e., compression stresses are positive, and snow elasticity is neglected.

5.1 Oedometric compression test

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 tests to
perform on snow. Under the lateral constraint E<sub>rr</sub> = 0, we
have

$$
\dot{E} = \begin{pmatrix} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & \dot{E}_{zz} \end{pmatrix}. \qquad (36)
$$

The quasi-static equilibrium div Σ = 0 implies that Σθθ = Σrr. Consequently, the macroscopic stress tensor is written

$$
\Sigma = \begin{pmatrix} \Sigma_{rr} & 0 & 0 \\ 0 & \Sigma_{rr} & 0 \\ 0 & 0 & \Sigma_{zz} \end{pmatrix}. \qquad (37)
$$

As a result, the two first strain rates and stress invariants are
written

$$
\begin{align}
E_1 &= \dot{E}_{zz}, \quad \bar{E}_2 = \sqrt{\frac{2}{3}} \dot{E}_{zz}, & S_1 &= 2\Sigma_{rr} + \Sigma_{zz}, \nonumber \\
\bar{S}_2 &= \sqrt{\frac{2}{3}}(\Sigma_{zz} - \Sigma_{rr}), & \text{with } \Sigma_{zz} &\ge \Sigma_{rr}. \tag{38}
\end{align}
$$

In this particular case, from Eqs. (33) and (34), it can be
shown that the lateral constraint $E_{rr} = 0$ implies that

$$
\frac{\Sigma_{rr}}{\Sigma_{zz}} = \frac{c(\phi) - 2f(\phi)}{c(\phi) + 4f(\phi)}, \qquad (39)
$$

<figure><img src="image_20.png" /><figcaption>**Figure 9.** Oedometric compression test (a) and triaxial compression test (b).</figcaption></figure>

and consequently

$$
\begin{align}
\Sigma_{\text{eq}} &= \sqrt{\frac{9c(\phi)f(\phi)}{4f(\phi) + c(\phi)}} \Sigma_{zz} \quad \text{and} \nonumber \\
\dot{E}_{zz} &= A \left( \frac{9f(\phi)c(\phi)}{4f(\phi) + c(\phi)} \right)^{\frac{n+1}{2}} \Sigma_{zz}^n. \tag{40}
\end{align}
$$

In oedometric experimental tests, the lateral pressure Σrr
is not easily accessible and it is often tempting to neglect
this pressure and interpret any oedometric compression test
as a uniaxial compression test. The relation (Eq. 39) can be
used to assess the relative importance of the confining pres-
sure with respect to the vertical stress. The evolution of this
ratio with respect to snow density is shown in Fig. 10 for
n ∈ {2, 3, 4.5}. It should be noticed that this ratio does not
depend on the axial strain rate ˙Ezz. As expected, this ratio in-
creases with increasing snow density and tends towards one
for φ = 0 due to the incompressibility of the ice skeleton.
In the whole range of snow compacities under consideration
/materialized by the gray zone in Fig. 10) the lateral pressure
represents 30 to 50 % of the vertical stress and cannot be ne-
glected in practice. Finally, let us remark that the evolution of
this ratio is similar to the one measured by Geindreau et al.
(1999a) and Viot and Stutz (2002) on metallic powders. Fig-
ure 10 also shows that this ratio is almost independent of the
exponent n, which is consistent with the experimental results
of Viot and Stutz (2002).

5.2 Triaxial compression test

During a triaxial test, a cylindric snow sample is submitted
simultaneously to an axial stress Σ<sub>zz</sub> and a lateral confin-
ing pressure Σ<sub>rr</sub>. The static equilibrium div Σ = 0 implies
that Σ<sub>θθ</sub> = Σ<sub>rr</sub>; consequently the macroscopic stress tensor
is written

$$
\Sigma = \begin{pmatrix} \Sigma_{rr} & 0 & 0 \\ 0 & \Sigma_{rr} & 0 \\ 0 & 0 & \Sigma_{zz} \end{pmatrix}. \qquad (41)
$$
<figure><img src="image_21.png" /><figcaption>**Figure 10.** Evolution of the ratio between the lateral and axial stresses with respect to snow relative density during any oedometric compression test for different values of the exponent *n* (numerical points as diamond points; functions corresponding to Eq. 29 as solid lines). The relative density range for typical snow samples is materialized by the gray zone.</figcaption></figure>

In these conditions, the constitutive Eq. (33) implies that $\dot{E}_{\theta\theta} = \dot{E}_{rr}$, and thus

$$ \dot{E} = \begin{pmatrix} \dot{E}_{rr} & 0 & 0 \\ 0 & \dot{E}_{rr} & 0 \\ 0 & 0 & \dot{E}_{zz} \end{pmatrix}. \qquad (42) $$

As a result, the two first strain rate and stress invariants are written

$$ \begin{aligned} E_1 &= 2\dot{E}_{rr} + \dot{E}_{zz}, & \bar{E}_2 &= \sqrt{\frac{2}{3}}(\dot{E}_{zz} - \dot{E}_{rr}), \\ S_1 &= 2\Sigma_{rr} + \Sigma_{zz}, & \bar{S}_2 &= \sqrt{\frac{2}{3}}|\Sigma_{zz} - \Sigma_{rr}|. \end{aligned} \qquad (43) $$

In this particular case, from Eqs. (33) and (34) it can be shown that

$$ \Sigma_{eq} = \sqrt{f(\phi)(2\Sigma_{rr} + \Sigma_{zz})^2 + c(\phi)(\Sigma_{zz} - \Sigma_{rr})^2} \qquad (44) $$

and

$$ \left\{ \begin{aligned} \dot{E}_{rr} &= A \Sigma_{eq}^{n-1} \left[ \left( 2f(\phi) + \frac{1}{2}c(\phi) \right) \Sigma_{rr} \right. \\ & \left. \qquad + \left( f(\phi) - \frac{1}{2}c(\phi) \right) \Sigma_{zz} \right] , \\ \dot{E}_{zz} &= A \Sigma_{eq}^{n-1} \left[ \left( 2f(\phi) - c(\phi) \right) \Sigma_{rr} \right. \\ & \left. \qquad + \left( f(\phi) + c(\phi) \right) \Sigma_{zz} \right] \end{aligned} \right. \qquad (45) $$

In the case of a uniaxial compression test, $\Sigma_{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), Fig. 11 presents the evolution of the snow densification rate (for $n = 4.5$ and $A = 1.5 \times 10^{-3} \text{ MPa}^{-n} \text{ s}^{-1}$) given by $E_1 = \dot{\rho}_{\text{snow}} / \rho_{\text{snow}}$ with respect to snow relative density when constant stresses are applied on the sample. As expected, this figure shows the following:

- Whatever the loading, the densification rate strongly decreases with increasing snow density. In the investigated range, i.e., $\rho_{\text{snow}} / \rho_{\text{ice}} \in [0.1, 0.6]$, the densification rate decreases by 9 orders of magnitude from $10^{-1}$ to $10^{-10} \text{ s}^{-1}$.

- For a given snow density, the loading conditions influence strongly the densification rate. Typically, when $\Sigma_{zz} = 10 \text{ kPa}$ the densification rate decreases by nearly 1 order of magnitude if the confining pressure $\Sigma_{rr}$ is reduced from 10 kPa (isotropic compression) to 0 kPa (uniaxial compression). In contrast, when $\Sigma_{rr} = 10 \text{ kPa}$, the densification rate increases by nearly 1 order of magnitude if the axial stress $\Sigma_{zz}$ increases from 10 kPa (isotropic compression) to 20 kPa (triaxial compression). As expected, this last result shows the increase in the densification rate with the increase in the
<figure><img src="image_22.png" /><figcaption>**Figure 11.** Predicted densification rate (for $n = 4.5$ and $A = 1.5 \times 10^{-3} \text{ MPa}^{-n} \text{ s}^{-1}$) with respect to snow relative density for four classical laboratory tests: an oedometric compression test (dash-dotted line), a uniaxial compression test (solid line), a triaxial compression test (dashed line) and an isotropic compression test (dotted line). The inset plot provides a zoom on the classical range of snow densities observed experimentally (gray background).</figcaption></figure>

deviatoric stress (i.e., $\bar{S}_2$). Even if the lateral confining pressure cannot be neglected during the oedometric test as highlighted in Fig. 10, the oedometric compression test results in a similar densification rate as the uniaxial compression test for the same axial stress $\Sigma_{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 dilating, which is not the case for the uniaxial compression test. Overall the two effects cancel out each other. Finally, above the classical snow relative density range ($\rho_{\text{snow}}/\rho_{\text{ice}} \ge 0.6$), the densification rate dramatically decreases for the oedometric and isotropic compression tests due to the ice incompressibility. As already underlined in Fig. 10 for 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. Figure 12 presents the evolution of the stress $\Sigma_{zz}$ versus snow relative density $\rho_{\text{snow}}/\rho_{\text{ice}} = 1 - \phi$ for two different values of strain rates $\dot{E}_{zz} \in \{10^{-7}; 10^{-5}\} \text{ s}^{-1}$ and the different loading

conditions (uniaxial, oedometric, isotropic and triaxial tests). This figure suggests the following:

- As expected, for a given strain rate, the stress $\Sigma_{zz}$ increases with increasing snow density.

- For a given strain rate and a given density, the stress $\Sigma_{zz}$ increases with increasing the lateral pressure $\Sigma_{rr}$ around the sample.

- For a given snow density, the stress $\Sigma_{zz}$ strongly increases with increasing strain rate, which is in accordance with the power law relationship.

- The ice viscoplastic behavior is recovered when $\rho_{\text{snow}}/\rho_{\text{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 ice incompressibility, the axial oedometer stress $\Sigma_{zz}$ tends towards $+\infty$.

In order to quantitatively compare the predictions of our model against the experimental results of Bartelt and
<figure><img src="image_23.png" /><figcaption>**Figure 12.** Predicted axial stress level (for $n = 4.5$ and $A = 1.5 \times 10^{-3}$ MPa$^{-n}$ s$^{-1}$) with respect to snow relative density for two imposed strain rates (line thickness) and for four classical laboratory tests: an oedometric compression test (dash-dotted line), a uniaxial compression test (solid line), a triaxial compression test (dashed line) and an isotropic compression test (dotted line). The inset plot provides a zoom on the classical range of snow densities observed experimentally (gray background).</figcaption></figure>

von Moos (2000), we consider a snow of density $\rho_{snow} = 255$ kg m$^{-3}$ (corresponding to $\rho_{snow}/\rho_{ice} = 0.27$) subjected to a confining pressure $\Sigma_{rr} = 2.5$ kPa and a strain rate $\dot{E}_{zz} = 2.2 \times 10^{-5}$ s$^{-1}$. In this case, the axial stress predicted by our model is $\Sigma_{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 a uniaxial compression test ($\Sigma_{rr} = 0$). In this case, the axial stress simply reads

$$ \Sigma_{zz}(\dot{E}_{zz}, \phi) = \left( \frac{\dot{E}_{zz}}{A(f(\phi) + c(\phi))^{\frac{n+1}{2}}} \right)^{\frac{1}{n}}. \quad (46) $$

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:

$$ \begin{aligned} \alpha_{\eta}(\phi) &= \frac{\Sigma_{zz}^{\text{snow}}}{(1-\phi)\Sigma_{zz}^{\text{ice}}} = \frac{\Sigma_{zz}(\dot{E}_{zz}, \phi)}{(1-\phi)\Sigma_{zz}(\dot{E}_{zz}, 0)} \\ &= \frac{1}{1-\phi} \left( \frac{1}{f(\phi) + c(\phi)} \right)^{\frac{n+1}{2n}}. \end{aligned} \quad (47) $$

This parameter compares the axial stress that a given snow sample can transmit ($\Sigma_{zz}^{\text{snow}}$) to a rough estimate of this stress given as a fraction of the axial stress transmitted in the case of ice ($(1-\phi)\Sigma_{zz}^{\text{ice}}$). In Fig. 13, the above theoretical expression of $\alpha_{\eta}(\phi)$ is compared with the experimental fit $\alpha_{\eta} = 0.0028 \exp(0.008\rho_{ice}(1-\phi))$ proposed by Bartelt and von Moos (2000). As expected, $\alpha_{\eta}(\phi)$ increases with increasing snow density. By definition, $\alpha_{\eta}(\phi)$ should vary between 0 and 1. We can observe that the theoretical expression of $\alpha_{\eta}(\phi)$ is strictly greater than 1 for $\rho_{snow}/\rho_{ice} \in [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 Sect. 4.2. An implicit relation between these parameters could help in order to en-
<figure><img src="image_24.png" /><figcaption>**Figure 13.** Comparison between our theoretical prediction for $\alpha_{\eta}$ (Eq. 47) and the experimental fit proposed by Bartelt and von Moos (2000). The points computed directly from the values for *f* and *c* reported in Table 3 are shown as diamond points. The relative density range for snow samples under consideration is materialized by the gray zone.</figcaption></figure>

sure that $\alpha_{\eta}(\phi)$ remains lower than 1 in the whole relative density range. Nevertheless, in the range of snow densities under consideration ($\rho_{\text{snow}}/\rho_{\text{ice}} \in [0.1, 0.6]$), $\alpha_{\eta}$ increases monotonously between roughly 0.1 and 0.6. This prediction is higher than the experimental fit proposed by Bartelt and von Moos (2000) (see Fig. 13). However, during this experiment, the macroscopic mechanical response probably results from both viscous deformation of ice skeleton and ruptures of ice bridges between snow grains. Even if our model is able to account for some microstructure modifications through the porosity dependence of the parameters *f* and *c*, the changes induced by the experimental testing conditions might exceed the scope of application of our model.

## 6 Conclusions

Despite the nonlinearity of the viscous ice constitutive equation, the image-based homogenization approach introduced by Wautier et al. (2015) was successfully adapted to the numerical homogenization of snow viscous behavior. It allows the viscoplastic response of any snow sample being computed from its X-ray tomographic image. In contrast to the elastic case, the macroscopic stress response is not a linear function of the imposed macroscopic strain anymore. As a result,

snow macroscopic response was investigated in terms of iso-dissipation curves in the planes of the two first strain rate and stress invariants. The shape and size of these curves characterize the strong coupling between the snow microstructure and the viscous ice behavior at the microscopic scale (power law with an exponent *n*). Different values for *n* were considered in this study to give some insight on their influence on the homogenized viscoplastic behavior of snow. Thanks to a few selected loading paths, an Abouaf model was fitted onto numerical simulation results. This formulation seems to be relevant to describe the snow viscoplastic behavior in the whole range of snow porosity under consideration ($0.43 < \phi < 0.87$), provided that snow microstructure is isotropic. For a given value of *n*, the influence of the snow microstructure on its viscoplastic response is described at the first order through two material functions *f*($\phi$) and *c*($\phi$) (see Eq. 29) depending on the porosity only.

The robustness of this Abouaf formulation was 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 and contrary to the case of ice, the ability of snow to exhibit a viscous behavior even under isotropic strain loading is recovered. The scope of 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 recover the experimental results of Desrues et al. (1980), Bartelt and von Moos (2000), Moos et al. (2003) and 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, 2016; Wautier et al., 2015) but the extension of our homogenized viscoplastic formulation to anisotropic snow types is quite challenging as the dimension of the invariant space will increase dramatically (Boehler, 1978; Liu, 1982; Hansen et al., 1991).

Finally, in the present work based on FEM simulations, the ice skeleton is viewed as a continuous polycrystalline material. The proposed methodology to identify and formulate the 3-D viscoplastic behavior of snow could be also applied to DEM simulations. For that purpose, the identification of the shape and crystalline orientations of every ice grain (Roland du Roscoat et al., 2011) as well as the knowledge of the viscoplastic contact laws (Burr et al., 2015b, a) would be of primary interest.

*Data availability.* The original 3-D images used in this study are further described in the Supplement of Calonne et al. (2012) and Wautier et al. (2015). They can be obtained on request (fREDERIC.FLIN@METEO.FR) after signing an agreement with Météo-France. More details on the numerical tools used to convert binary images into meshes can be found in Fang and Boas (2009) and they can be downloaded at http://iso2mesh.sourceforge.net/cgi-bin/index.cgi. More details on the plug-in used to perform the homogenization test (KUBC) can be found in Lejeunes and Bourgeois (2011) and they can be downloaded at http://www.lma.cnrs-mrs.fr/spip.php?article171.

*Competing interests.* The authors declare that they have no conflict of interest.

*Acknowledgements.* We thank Stéphane Lejeunes and Stéphane Bourgeois from the LMA for sharing with us the code of their plug-in Homtools and for helping us to use it in our scripts. Tomographic images have been acquired at ESRF (ID19) and 3SR Lab. CNRM/CEN is part of the LabEx OSUG@2020 (ANR10 LABX56). 3SR lab is part of the LabEx Tec 21 (Investissements d'Avenir – grant agreement ANR11 LABX0030).

The authors gratefully acknowledge the anonymous reviewers, Andrew Hansen and Maurine Montagnat for their comments and advice, which helped improve the quality of the final paper.

Edited by: Martin Schneebeli

Reviewed by: Andrew C. Hansen and two anonymous referees

## References

Abouaf, M.: Modélisation de la compaction de poudres métalliques frittées, PhD thesis, USMG-INP, Grenoble, France, 1985.

Abouaf, M. and Chenot, J.: Finite element simulation of hot isostatic pressing of metal powders, Int. J. Numer. Meth. Eng., 25, 191–212, 1988.

Adams, E. E. and Walters, D. J.: Fine structure layering in radiation recrystallized snow, in: International Snow Science Workshop 2014 Proceedings, 29 September–3 October 2014, Banff, Canada, 29–34, ISSW, 2014.

Auriault, J.-L., Bouvard, D., Dellis, C., and Lafer, M.: Modeling of hot compaction of metal powder by homogenization, Mech. Mater., 13, 247–275, 1992.

Auriault, J.-L., Royer, P., and Geindreau, C.: Filtration law for power law fluids in anisotropic media, Int. J. Eng. Sci., 40, 1151–1163, 2002.

Auriault, J.-L., Boutin, C., and Geindreau, C.: Homogenization of coupled phenomena in heterogenous media, vol. 149, John Wiley & Sons, London, UK, 2010.

Bartelt, P. and von Moos, M.: Triaxial tests to determine a microstructure-based snow viscosity law, Ann. Glaciol., 31, 457–462, 2000.

Boehler, J.: Lois de comportement anisotrope des milieux continus, J. Mécanique, 17, 153–190, 1978.

Brzoska, J.-B., Coléou, C., Lesaffre, B., Borel, S., Brissaud, O., Ludwig, W., Boller, E., and Baruchel, J.: 3D visualization of snow samples by microtomography at low temperature, ESRF Newsletter, 32, 22–23, 1999.

Burr, A., Philip, A., and Martin, C.: Etude expérimentale de la déformation viscoplastique de cylindres monocristallins de glace, in: Plasticité 2015, 28–30 April 2015, Autrans, France, 2015a.

Burr, A., Trecourt, P., Philip, A., and Martin, C.: Densification of firn using the Discrete Element Method, in: ESMC 2015, 9th European Solid Mechanics Conference, 6–10 July 2015, Madrid, Spain, 2015b.

Burr, A., Noël, W., Trecourt, P., Bourcier, M., Gillet-Chaulet, F., Philip, A., and Martin, C. L.: The anisotropic contact response of viscoplastic monocrystalline ice particles, Acta Mater., 132, 576–585, 2017.

Calonne, N., Geindreau, C., Flin, F., Morin, S., Lesaffre, B., Roland du Roscoat, S., and Charrier, P.: 3-D image-based numerical computations of snow permeability: links to specific surface area, density, and microstructural anisotropy, The Cryosphere, 6, 939–951, https://doi.org/10.5194/tc-6-939-2012, 2012.

Calonne, N., Flin, F., Geindreau, C., Lesaffre, B., and Rolland du Roscoat, S.: Study of a temperature gradient metamorphism of snow from 3-D images: time evolution of microstructures, physical properties and their associated anisotropy, The Cryosphere, 8, 2255–2274, https://doi.org/10.5194/tc-8-2255-2014, 2014.
Calonne, N., Flin, F., Lesaffre, B., Dufour, A., Roulle, J., Puglièse, P., Philip, A., Lahoucine, F., Geindreau, C., Panel, J.-M., Rolland du Roscoat, S., and Charrier, P.: CellDyM: A room temperature operating cryogenic cell for the dynamic monitoring of snow metamorphism by time-lapse X-ray microtomography, Geophys. Res. Lett., 42, 3911–3918, https://doi.org/10.1002/2015GL063541, 2015.

Chandel, C., Srivastava, P. K., and Mahajan, P.: Micromechanical analysis of deformation of snow using X-ray tomography, Cold Reg. Sci. Technol., 101, 14–23, 2014.

Chen, S. and Baker, I.: Evolution of individual snowflakes during metamorphism, J. Geophys. Res., 115, D21114, https://doi.org/10.1029/2010JD014132, 2010.

Cresseri, S. and Jommi, C.: Snow as elastic viscoplastic bonded continuum: a modelling approach, Rivista Italiana Di Geotecnica, 4, 43–58, 2005.

Cresseri, S., Genna, F., and Jommi, C.: Numerical integration of an elastic-viscoplastic constitutive model for dry metamorphosed snow, Int. J. Numer. Anal. Met., 34, 1271–1296, 2009.

Danas, K., Idiart, M., and Ponte Casteneda, P.: A homogenization-based constitutive model for isotropic viscoplastic porous media, Int. J. Solids Struct., 45, 3392–3409, 2008.

Desrues, J., Darve, F., Flavigny, E., Navarre, J., and Taillefer, A.: An incremental formulation of constitutive equations for deposited snow, J. Glaciol., 25, 289–307, 1980.

Dormieux, L. and Bourgeois, E.: Introduction à la micromécanique des milieux poreux, Presses de l'École nationale des ponts et chaussées, Paris, France, 2002.

Duva, J. and Crow, P.: The densification of powders by power-law creep during hot isostatic pressing, Acta Metall. Mater., 40, 31–35, 1992.

Fang, Q. and Boas, D. A.: Tetrahedral mesh generation from volumetric binary and grayscale images, in: ISBI'09. IEEE International Symposium on Biomedical Imaging: From Nano to Macro, 2009, 28 June–1 July 2009, Melbourne, Australia, 1142–1145, IEEE, available at: http://iso2mesh.sourceforge.net/cgi-bin/index.cgi (last access: 20 June 2017), 2009.

Fierz, C., Armstrong, R. L., Durand, Y., Etchevers, P., Greene, E., McClung, D. M., Nishimura, K., Satyawali, P. K., and Sokratov, S. A.: The International Classification for Seasonal Snow on the Ground, UNESCO/IHP, Paris, France, 2009.

Flin, F. and Brzoska, J.-B.: The temperature-gradient metamorphism of snow: vapour diffusion model and application to tomographic images, Ann. Glaciol., 49, 17–21, 2008.

Fritzen, F., Forest, S., Bohlke, T., Kondo, D., and Kanit, T.: Computational homogenization of elasto-plastic porous metals, Int. J. Plasticity, 29, 102–119, 2012.

Geindreau, C. and Auriault, J.-L.: Investigation of the mechanical behaviour of alloys in the semi-solid state by homogenization, Mech. Mater., 31, 535–551, 1999.

Geindreau, C., Bouvard, D., and Doremus, P.: Constitutive behaviour of metal powder during hot forming. Part I: experimental investigation with lead powder as a simulation material, Eur. J. Mech. A-Solids, 18, 581–596, 1999a.

Geindreau, C., Bouvard, D., and Doremus, P.: Constitutive behaviour of metal powder during hot forming. Part II: Unified viscoplastic modelling, Eur. J. Mech. A-Solids, 18, 597–615, 1999b.

Green, R.: A plasticity theory for porous solids, Int. J. Mech. Sci., 14, 215–224, 1972.

Hagenmuller, P., Theile, T. C., and Schneebeli, M.: Numerical simulation of microstructural damage and tensile strength of snow, Geophys. Res. Lett., 41, 86–89, https://doi.org/10.1002/2013GL058078, 2014.

Hagenmuller, P., Chambon, G., and Naaim, M.: Microstructure-based modeling of snow mechanics: a discrete element approach, The Cryosphere, 9, 1969–1982, https://doi.org/10.5194/tc-9-1969-2015, 2015.

Hansen, A. C., Blackketter, D. M., and Walrath, D. E.: An invariant-based flow rule for anisotropic plasticity applied to composite materials, J. Appl. Mech., 58, 881–888, https://doi.org/10.1115/1.2897701, 1991.

Johnson, J. B. and Hopkins, M. A.: Identifying microstructural deformation mechanisms in snow using discrete element modeling, J. Glaciol., 51, 432–442, 2005.

Kaempfer, T. U., Schneebeli, M., and Sokratov, S.: A microstructural approach to model heat transfer in snow, Geophys. Res. Lett., 32, L21503, https://doi.org/10.1029/2005GL023873, 2005.

Kanit, T., Forest, S., Galliet, I., Mounoury, V., and Jeulin, D.: Determination of the size of the representative volume element for random composites: statistical and numerical approach, Int. J. Solids Struct., 40, 3647–3679, 2003.

Köchle, B. and Schneebeli, M.: Three-dimensional microstructure and numerical calculation of elastic properties of alpine snow with a focus on weak layers, J. Glaciol., 60, 705–713, 2014.

Lejeunes, S. and Bourgeois, S.: Une Toolbox Abaqus pour le calcul de propriétés effectives de milieux hétérogènes, in: 10e colloque national en calcul des structures (CSMA 2011), 9–13 May 2011, Presqu'île de Giens, France, available at: http://www.lma.cnrs-mrs.fr/spip.php?article171 (last access: 20 June 2017), 2011.

Lemaitre, J. and Chaboche, J.-L.: Mécanique des Matériaux Solides, Dunod, Paris, France, 1985.

Liu, I.-S.: On representations of anisotropic invariants, International Journal of Engineering Science, 20, 1099–1109, 1982.

Löwe, H., Riche, F., and Schneebeli, M.: A general treatment of snow microstructure exemplified by an improved relation for thermal conductivity, The Cryosphere, 7, 1473–1480, https://doi.org/10.5194/tc-7-1473-2013, 2013.

Mede, T., Chambon, G., Hagenmuller, P., and Nicot, F.: Mechanical behaviour of weak snow layers: modelling a porous structure of sintered grains, Powder and grains 2017 – 8th International Conference on Micromechanics of Granular Media, 3–7 July 2017, Montpellier, France, https://doi.org/10.1051/epjconf/201714006008, 2017.

Mellor, M.: A review of basic snow mechanics, US Army Cold Regions Research and Engineering Laboratory, Snow mechanics Symposium, April 1974, Grindelwald, Switzerland, 1974.

Moos, M. V., Bartelt, P. A., Zweidler, A., and Bleiker, E.: Triaxial tests on snow at low strain rate. Part I: Experimental device, J. Glaciol., 49, 81–90, 2003.

Narita, H.: An experimental study on tensile fracture of snow, Contributions from the institute of Low Temperature Science, 32, 1–37, 1984.

Navarre, J. P., Meyssonnier, J., and Vagnon, A.: 3D numerical model of snow deformation without failure and its application to cold room mechanical tests, Cold Reg. Sci. Technol., 50, 3–12, 2007.
Orgéas, L., Geindreau, C., Auriault, J.-L., and Bloch, J.-F.: Upscal-
ing the flow of generalised Newtonian fluids through anisotropic
porous media, J. Non-Newton. Fluid, 145, 15–29, 2007.

Pieritz, R. A., Brzoska, J.-B., Flin, F., Lesaffre, B., and Coléou,
C.: From snow X-ray microtomograph raw volume data to mi-
cromechanics modeling: first results, Ann. Glaciol., 38, 52–58,
https://doi.org/10.3189/172756404781815176, 2004.

Pinzer, B. R., Schneebeli, M., and Kaempfer, T. U.: Va-
por flux and recrystallization during dry snow metamor-
phism under a steady temperature gradient as observed by
time-lapse micro-tomography, The Cryosphere, 6, 1141–1155,
https://doi.org/10.5194/tc-6-1141-2012, 2012.

Rolland du Roscoat, S., King, S., Philip, A., Reischig, A., Ludwig, P., Flin, F., and Meysonnier, J.: Analysis of snow microstructure by means of X-ray diffraction contrast tomography, Adv. Eng. Mater., 13, 128–135, 2011.

Salm, B.: Mechanical properties of snow, Rev. Geophys., 20, 1–19,
1982.

Sanchez, L., Ouedraogo, E., Federzoni, L., and Stutz, P.: New vis-
coplastic model to simulate hot isostatic pressing, Powder Met-
all., 45, 329–334, 2002.

Scapozza, C. and Bartelt, P.: Triaxial tests on snow at low strain rate. Part II. Constitutive behaviour, J. Glaciol., 49, 91–101, 2003.

Schleef, S., Löwe, H., and Schneebeli, M.: Hot-pressure sintering of low-density snow analyzed by X-ray microtomography and in situ microcompression, Acta Mater., 71, 185–194, 2014.

Schneebeli, M.: Numerical simulation of elastic stress in the mi-
crostructure of snow, Ann. Glaciol., 38, 339–342, 2004.

Schulson, E. M. and Duval, P.: Creep and fracture of ice, Cambridge
University Press, Cambridge, UK, 2009.

Schweizer, J. and Camponovo, C.: The temperature dependence of the effective elastic shear modulus of snow, Cold Reg. Sci. Technol., 35, 55–64, 2002.

Scopozza, C. and Bartelt, P. A.: The influence of temperature on the small strain viscous deformation mechanics of snow: comparison with polycrystalline ice, Ann. Glaciol., 37, 90–96, 2003a.

Scopozza, C. and Bartelt, P. A.: Triaxial tests on snow at low strain rate. Part II: Constitutive modelling, J. Glaciol., 49, 91–101, 2003b.

Shapiro, L. H., Johnson, J. B., Sturm, M., and Blaisdell, G. L.: Snow mechanics: review of the state of knowledge and applications ap- plications (No. CRREL-97-3), Cold Regions Research And En- gineering Lab, Hanover, NH, USA, 1997.

Smith, G.: On Isotropic Functions of Symmetric Tensors Skew-Symmetric Tensors and Vectors, Int. J. Eng. Sci., 9, 899–916, 1971.

Sofronis, P. and McMeeking, R.: Creep of power-law material con-
taining spherical voids, J. Appl. Mech., 59, S88–S95, 1992.

Srivastava, P., Mahajan, P., Satyawali, P., and Kumar, V.: Obser-
vation of temperature gradient metamorphism in snow by X-
ray computed microtomography: measurement of microstruc-
ture parameters and simulation of linear elastic properties, Ann.
Glaciol., 51, 73–82, 2010.

Srivastava, P. K., Chandel, C., Mahajan, P., and Pankaj, P.: Predic-
tion of anisotropic elastic properties of snow from its microstruc-
ture, Cold Reg. Sci. Technol., 125, 85–100, 2016.

Storakers, B., Fleck, N. A., and McMeeking, R. M.: The viscoplas-
tic compaction of composite powders, J. Mech. Phys. Solids, 47,
785–815, 1999.

Suquet, P.: Overall potentials and extremal surfaces of power law or ideally plastic composites, J. Mech. Phys. Solids, 41, 981–1002, 1993.

Theile, T., Lowe, H., Theile, T. C., and Schneebeli, M.: Simulating creep of snow based on microstructure and anisotropic deformation of ice, Acta Mater., 59, 7104–7113, 2011.

Viot, P. and Stutz, P.: Nouveau dispositif expérimental pour l'étude du comportement viscoplastique des poudres métalliques à hautes températures: application à une poudre de cuivre, C. R. Mecanique, 330, 653–659, 2002.

Wang, X. and Baker, I.: Observation of the microstructural evolu-
tion of snow under uniaxial compression using X-ray computed
microtomography, J. Geophys. Res.-Atmos., 118, 12–371, 2013.

Wautier, A., Geindreau, C., and Flin, F.: Linking snow microstruc-
ture to its macroscopic elastic stiffness tensor: A numeri-
cal homogenization method and its application to 3-D images
from X-ray tomography, Geophys. Res. Lett., 42, 8031–8041,
https://doi.org/10.1002/2015GL065227, 2015.