A model for French-press experiments of dry snow compaction

Snow densification stores water in alpine regions and transforms snow into ice on the surface of glaciers. Despite its importance in determining snow-water equivalent and glacier-induced sea level rise, we still lack a complete understanding of the physical mechanisms underlying snow compaction. In essence, compaction is a rheological process, where the rheology evolves with depth due to variation in temperature, pressure, humidity, and meltwater. The rheology of snow compaction can be determined in a few ways, for example, through empirical investigations (e.g., Herron and Langway, 1980), by microstructural considerations (e.g., Alley, 1987), or by measuring the rheology directly, which is the approach we take here. Using a “French-press” or “cafetière-à-piston” compression stage, Wang and Baker (2013) compressed numerous snow samples of different densities. Here we derive a mixture theory for compaction and airflow through the porous snow to compare against these experimental data. We find that a plastic compaction law explains experimental results. Taking standard forms for the permeability and effective pressure as functions of the porosity, we show that this compaction mode persists for a range of densities and overburden loads. These findings suggest that measuring compaction in the lab is a promising direction for determining the rheology of snow through its many stages of densification.

1 Introduction

Snow densification in alpine and polar regions transforms snowflakes into ice crystals. On the surface of glaciers and ice sheets, fresh snow is buried by new snow each winter, thereby slowly transforming into firn and then glacial ice as it compresses and descends. In cold and dry environments (e.g., melt-free accumulation areas of mountain glaciers, interior Greenland, and Antarctica), surface snow evolves due to temperature gradient metamorphism and atmospheric interactions (Dadic et al., 2008; Chen and Baker, 2010b). In the region below the surface, snow compaction is thought to occur in three stages (Herron and Langway, 1980; Arnaud et al., 2000; Cuffey and Paterson, 2010). In the first stage near the surface, compaction occurs by grain growth and rearrangement due to grain boundary sliding (Alley, 1987). In the second stage, compaction occurs through increasing overburden pressure inducing creep deformation and sintering (Wilkinson and Ashby, 1975; Wilkinson, 1988; Arnaud et al., 1998; Spencer et al., 2001). In the final stage, interconnected pores have closed off, and further compaction is caused by air bubble compression (Alley and Bentley, 1988; Salamatin et al., 1997; Gregory et al., 2014). In wet environments, such as the percolation zone of mountain glaciers and Greenland as well as many ice shelves of Antarctica, compaction occurs by a combination of dry snow compaction processes and refreezing of meltwater, which can either enhance or detract from the densification processes just mentioned (Colbeck, 1976; Machguth et al., 2016; Meyer and Hewitt, 2017). Although meltwater percolation is an important part of the compaction process in many areas (e.g., Colbeck, 1972; Bartelt and Lehning, 2002; Wever et al., 2014; Steger et al., 2017), we will only consider dry snow compaction here.

An important reason for studying snow compaction on glaciers, ice sheets, and snowfields is to connect a change in surface elevation to a volume of stored water. The total water volume stored in glaciers and snowpacks is important to know for current as well as future water resources and sea level rise considerations. Additionally, compaction is important for ice core analysis (Barnola et al., 1991; Goujon et al., 2003; Cuffey and Paterson, 2010). At the bottom of the firn column, the difference between the ice age and the gas age is an important input into paleoclimate reconstructions of temperature in ice cores and must be estimated using a compaction model (Arnaud et al., 2000). Compaction can be measured through tracking relative displacements of metal markers (Hawley et al., 2004) or snow features in optical stratigraphy (Hawley and Waddington, 2011), autonomous phase-sensitive radio-echo sounding (Nicholls et al., 2015), vertical strain fiber-optic sensors (Zumberge et al., 2002), the relative-displacement “coffee can” technique (Hamilton et al., 1998; Hamilton and Whillans, 2000), and a continuously recording coffee can method (Arthern et al., 2010; Stevens, 2018).

In one-dimensional compaction, the body force acting on the snow is gravity, leading to an overburden pressure σ_o, which increases with depth z below the surface z_s. Given mass and momentum conservation, all that is needed to predict the densification rate (i.e., surface ice velocity w_i(z_s) and strain rate within the snow column $\partial w_{\text{i}}/\partial z$) is a constitutive relationship between stress σ_o and strain rate $\partial w_{\text{i}}/\partial z$ in one vertical dimension. Thus, snow compaction can be thought of as a problem of rheology, where the rheology is complicated due to variation with depth, temperature, humidity, and water content, among many other physical processes. Three common approaches to characterizing the rheology of a material (Tanner, 2000) are empirical analysis (e.g., Herron and Langway, 1980), microstructural analysis (e.g., Alley, 1987), and experimental analysis. It is this last approach that we describe in this paper.

The standard compaction law applied to most glaciers and ice sheets is a one-dimensional relationship for the rate of change of snow density ρ with depth; i.e.,

$$\begin{array}{}\text{(1)}& {\displaystyle \frac{\partial \mathit{\rho }}{\partial t}}+w_{\text{i}}{\displaystyle \frac{\partial \mathit{\rho }}{\partial z}}=-\mathcal{C}\left(\mathit{\rho },{\mathit{\rho }}_{\text{i}},T,\dot{a},{\mathit{\sigma }}_{\mathrm{o}},{\displaystyle \frac{\partial w_{\text{i}}}{\partial z}},\mathrm{\dots }\right),\end{array}$$

where 𝒞 is the compaction function and can depend on the snow density ρ, ice density ρ_i, temperature T, accumulation rate $\dot{a}$, overburden pressure σ_o, vertical strain rate $\partial w_{\text{i}}/\partial z$, humidity, water saturation, grain size, and potentially many other physical processes (Spiegelman, 1993; Arthern et al., 2010; Lundin et al., 2017). Herron and Langway (1980) take $\mathcal{C}=c(\dot{a},T)\mathit{\rho }$, where c is an empirical function. However, including the accumulation rate in the compaction function is dubious as it really should enter the mass conservation equations as a boundary condition (Meyer and Hewitt, 2017). Other forms of the righthand side of Eq. (1) are discussed in Zwally and Li (2002), Reeh (2008), Morris and Wingham (2014), and Morris (2018). In this paper, we experimentally and mathematically analyze the function dependence of 𝒞 on overburden pressure and strain rate.

We now summarize the outline of the paper. In Sect. 2, we describe the laboratory compaction experiments of Wang and Baker (2013) and the samples used. In Sect. 3, we construct a continuum model to describe the experiments and implement our model numerically. Lastly, in Sect. 4, we compare the theoretical predictions for the snow compaction with the Wang and Baker (2013) data, showing that the theory does an excellent job explaining the snow compaction. We additionally discuss how this can be implemented into a compaction model such as Eq. (1) and our plans for future work, before a short conclusion section.

Table 1 Initial values of density and specific surface area for the four sintered low-temperature snow samples used in Wang and Baker (2013). Naming convention follows Fierz et al. (2009).

2 Experiments

To understand the microstructural origin of macroscopic snow material properties, the role of snow microstructure in avalanche initiation, and snow metamorphosis, Wang and Baker (2013) performed continuous compression tests on snow samples collected near Dartmouth College in Hanover, NH (USA). In their experiments, Wang and Baker (2013) focused on nine samples, ranging from (i) relatively warm, freshly fallen, low-density snow to (ii) snow that was collected during cold air temperatures (ca. −7 to −9 ^∘C) and placed in a −30 ^∘C cold room for 1 year, allowing for the snow crystals to sinter. In the sintering process, bonds form between snow crystals (Male, 1980; Chen and Baker, 2010a; Wang and Baker, 2017) and the snow evolved into round and well-connected snow grains during the year in the cold room (Chen and Baker, 2010b). In this paper, we focus on the four sintered samples from Wang and Baker (2013), as they are most applicable to alpine snowpacks and firn on glaciers. The naming convention for these sintered low-temperature (SLT) samples follows Fierz et al. (2009), where the highest density sample is SLT-1 and the lowest density sample is SLT-4 (see Table 1). Before compaction, the samples were vibrated to construct samples with different densities, undoubtedly breaking some of the sintered bounds. These sintered and vibrated samples represent a range of densities with relatively small differences in specific surface area (SSA; surface-to-volume ratio, viz. Table 1). The snow samples were then extruded as cylinders, 15.7 mm in diameter and 18 mm tall, and placed on a Skyscan material testing compression stage (i.e., a French-press or cafetière-à-piston-style coffee maker) inside of a microscopic computed tomography (micro CT) machine. A schematic of the compression stage is shown in Fig. 1. Both the top and bottom plates are impermeable but the air contained within the snow is able to flow out the side of the open sample. The samples were compressed at a constant rate of 12.7 mm h⁻¹ (i.e., 111 m yr⁻¹ or a pore-scale Reynolds number of $\mathit{Re}={\mathrm{10}}^{-\mathrm{4}}$ with 1 mm as a characteristic pore size) for the full 5.7 mm dynamic range of the Skyscan compression stage. To avoid termination effects, we only consider the first 5 mm of the range (Wang and Baker, 2013). The compression stage measured the load required to maintain a constant displacement rate for each snow sample. Wang and Baker (2013) then plot the loading stress as a function of the displacement, from which we can derive a stress-strain curve. In the next section, we develop a model for constant-displacement-rate compaction experiments, and in Sect. 4 we compare the model predictions for stress versus displacement against the Wang and Baker (2013) experimental data.

Figure 1 Schematic diagram of the snow French press: (a) idealized framework for the theory, including the boundary conditions; and (b) experimental setup from Wang and Baker (2013) showing the snow sample housed within the sample chamber and the upward motion of the lower plate. The sample is small enough that the effects of gravity can be ignored. The top and bottom plates are impermeable, but air can flow out of the sides of the snow sample.

3 Model

Compaction occurs in a variety of natural and industrial processes, such as sedimentary basin formation, paper pulp dewatering, and particle flocculation, and has motivated numerous experiments and mathematical models (e.g., Landman et al., 1991; Fowler and Noon, 1999; Fowler, 2011; Hewitt et al., 2016a, b; Paterson et al., 2019). In this section, we describe the theory outlined by Hewitt et al. (2016b), although we write the theory in terms of porosity ϕ rather than the solid fraction.

For dry snow that is composed of solely air and ice particles, snow density ρ is given as

$$\begin{array}{}\text{(2)}& \mathit{\rho }=\mathit{\varphi }{\mathit{\rho }}_{\text{a}}+(\mathrm{1}-\mathit{\varphi }){\mathit{\rho }}_{\text{i}},\end{array}$$

where ϕ is the porosity, i.e., the void space (Gray, 1996). The density of air is ρ_a and the density of pure ice is ρ_i, both of which we treat as constant. Although air density increases when pressurized in a confined space, in an open system where the air is able to flow, the air density is approximately constant. In the firn column of a glacier or ice sheet, the transition between freely flowing air and confined air is known as the pore-close-off depth. Below the pore-close-off depth the air density will increase with pressure, and above the pore-close-off depth, the air is approximately constant density.

Now from the expression in Eq. (2), it is clear that density variation in snow is equivalent to variation in porosity: near the surface of a glacier or snowpack, the snow air content is high and then density is closer to ρ_a, whereas at depth, the snow air content is low, and the density approaches ρ_i. Treating snow as a mixture of air and ice, mass conservation of each component may be written as

$$\begin{array}{}\text{(3)}& {\displaystyle }& {\displaystyle \frac{\partial \left({\mathit{\rho }}_{\text{a}}\mathit{\varphi }\right)}{\partial t}}+\mathbf{\nabla }\cdot \left({\mathit{\rho }}_{\text{a}}\mathit{\varphi }{\mathit{u}}_{\text{a}}\right)=m_{\text{s}},\text{(4)}& {\displaystyle }& {\displaystyle \frac{\partial \left[{\mathit{\rho }}_{\text{i}}(\mathrm{1}-\mathit{\varphi })\right]}{\partial t}}+\mathbf{\nabla }\cdot \left[{\mathit{\rho }}_{\text{i}}(\mathrm{1}-\mathit{\varphi }){\mathit{u}}_{\text{i}}\right]=-m_{\text{s}},\end{array}$$

where u_a is the air velocity, u_i is the ice velocity, and m_s is the mass exchange between each phase due to sublimation. Again, we treat the densities of each substance as constant, meaning that it is so easy for air to escape from the snow that the air density does not change. We also neglect phase change (i.e., m_s=0) and restrict our attention to one vertical dimension. Thus, we can write this model as

$$\begin{array}{}\text{(5)}& {\displaystyle }& {\displaystyle \frac{\partial \mathit{\varphi }}{\partial t}}+{\displaystyle \frac{\partial }{\partial z}}\left(\mathit{\varphi }{w}_a\right)=\mathrm{0},\text{(6)}& {\displaystyle }& {\displaystyle }-{\displaystyle \frac{\partial \mathit{\varphi }}{\partial t}}+{\displaystyle \frac{\partial }{\partial z}}\left[(\mathrm{1}-\mathit{\varphi })w_{\text{i}}\right]=\mathrm{0}.\end{array}$$

Now adding Eqs. (5) and (6) gives the insight that decreasing the porosity requires evacuating the incompressible air out of the pore space or, alternatively, air motion allows for snow compaction; i.e.,

$$\begin{array}{}\text{(7)}& {\displaystyle \frac{\partial }{\partial z}}\left[\mathit{\varphi }{w}_a+(\mathrm{1}-\mathit{\varphi })w_{\text{i}}\right]=\mathrm{0}.\end{array}$$

In other words, there is an exact volumetric trade-off where the pore space occupied by air is filled by ice during compaction.

In the Wang and Baker (2013) experiments, the sample has an initial height of $z=h_{\mathrm{0}}=\mathrm{18}$ mm (see Fig. 1), and, therefore, the overburden pressure is negligible when compared to the applied stress. We take the vertical coordinate system to increase downward. At the top of the sample z=0, the wall is impermeable; i.e., $w_a=w_{\text{i}}=\mathrm{0}$. Thus, we can integrate Eq. (7) to find that

$$\begin{array}{}\text{(8)}& w_{\text{i}}+\mathit{\varphi }\left(w_a-w_{\text{i}}\right)=\mathrm{0}.\end{array}$$

We model the airflow through porous snow using Darcy's law, which is given as

$$\begin{array}{}\text{(9)}& \mathit{\varphi }\left(w_a-w_{\text{i}}\right)=-{\displaystyle \frac{k\left(\mathit{\varphi }\right)}{\mathit{\mu }}}{\displaystyle \frac{\partial p}{\partial z}}\end{array}$$

for permeability k(ϕ), air viscosity μ, and air pressure p. As described by Hewitt et al. (2016b), Eq. (8) relates the solid ice velocity w_i to the Darcy velocity of the airflow through the pores. Thus, combining Eqs. (8) and (9), we find that

$$\begin{array}{}\text{(10)}& w_{\text{i}}={\displaystyle \frac{k\left(\mathit{\varphi }\right)}{\mathit{\mu }}}{\displaystyle \frac{\partial p}{\partial z}};\end{array}$$

that is, the ice particle velocity is determined by the pressure gradient driving air motion.

Force balance in the vertical direction is given as

$$\begin{array}{}\text{(11)}& {\displaystyle \frac{\partial \mathrm{\Sigma }}{\partial z}}=\mathrm{0},\end{array}$$

where Σ is the vertical bulk compressive stress (Hewitt et al., 2016b). We define effective pressure N as the difference between the compressive stress and air pressure, i.e., $N=\mathrm{\Sigma }-p$, and write Eq. (10) as

$$\begin{array}{}\text{(12)}& w_{\text{i}}=-{\displaystyle \frac{k\left(\mathit{\varphi }\right)}{\mathit{\mu }}}{\displaystyle \frac{\partial N}{\partial z}}.\end{array}$$

As is typical in soil mechanics (Tulaczyk et al., 2000) and other applications of compaction (Fowler, 2011), we can relate the effective pressure N to the porosity ϕ through an empirical constitutive relationship, so that N is identically equal to N(ϕ). This relationship is a form of plasticity as there is a one-to-one correspondence between load and density without reference to displacement, strain rate, or stress history. Inserting this into Eq. (12) and putting that into Eq. (6) gives

$$\begin{array}{}\text{(13)}& {\displaystyle \frac{\partial \mathit{\varphi }}{\partial t}}={\displaystyle \frac{\partial }{\partial z}}\left\{(\mathrm{1}-\mathit{\varphi })\left[-N^{\prime }\left(\mathit{\varphi }\right)\right]{\displaystyle \frac{k\left(\mathit{\varphi }\right)}{\mathit{\mu }}}{\displaystyle \frac{\partial \mathit{\varphi }}{\partial z}}\right\},\end{array}$$

which is a nonlinear diffusion equation for the porosity, and $N^{\prime }\left(\mathit{\varphi }\right)=\text{d}N/\text{d}\mathit{\varphi }$ is the derivative of the effective pressure with respect to the porosity.

Figure 2 Effective pressure N(ϕ) and permeability k(ϕ) constitutive relations as functions of the porosity ϕ for the snow compaction theory. The gray dots are permeability measurements of firn cores from Summit, Greenland (Adolph and Albert, 2014), and show excellent agreement with the Kozeny–Carman permeability model, Eq. (14) with a=3 and b=2.

3.1 Constitutive relations

Equation (13) is a general model for mechanical compaction, and the application to a specific system comes largely from the choice of constitutive relations for k(ϕ) and N(ϕ). A common choice for the permeability is the Kozeny–Carman model; i.e.,

$$\begin{array}{}\text{(14)}& k\left(\mathit{\varphi }\right)=k_{\mathrm{0}}{\displaystyle \frac{{\mathit{\varphi }}^a}{(\mathrm{1}-\mathit{\varphi }{)}^b}},\end{array}$$

which is commonly used throughout poromechanics (Rice and Cleary, 1976; McKenzie, 1984; Schoof and Hewitt, 2016) and has been extensively evaluated for snow and firn (Albert and Shultz, 2002; Adolph and Albert, 2013, 2014; Keegan et al., 2014), as shown in Fig. 2. Typical values for the exponents are a=3 and b=2. A competing model is a logarithmic permeability (Tait and Jaupart, 1992; Katz and Worster, 2008), of the form

$$\begin{array}{}\text{(15)}& k_{\mathrm{\ell }}\left(\mathit{\varphi }\right)=-k_{\mathrm{1}}{\mathit{\varphi }}^{\mathrm{2}}\ln \left(\mathrm{1}-\mathit{\varphi }\right),\end{array}$$

and there are many other permeability models (e.g., Hewitt et al., 2016b). Here we use the Kozeny–Carman model. However, before future compaction experiments, we will measure the permeability and determine the constitutive relation as well as associated parameters that best represent each sample.

The dependence of effective pressure on the porosity of a sample has not to our knowledge been measured for snow or firn. In other systems (e.g., Hewitt et al., 2016b), a plastic constitutive relationship between the effective pressure and porosity is given by

$$\begin{array}{}\text{(16)}& N\left(\mathit{\varphi }\right)=N_{\mathrm{0}}{\displaystyle \frac{(\mathrm{1}-\mathit{\varphi }{)}^n}{{\mathit{\varphi }}^m}},\end{array}$$

where typical exponents are n=3 and m=2 (see Fig. 2). An exponent of $n=\mathrm{3}/\mathrm{2}$ in the numerator is sometimes used for very porous materials (Hewitt et al., 2016b). Just as for the permeability, in the future we plan to measure the constitutive relation for effective pressure with porosity of each sample.

3.2 Boundary and initial conditions

As shown in Fig. 1, the samples were small and loaded from the bottom. The mathematical model for mechanical compaction results in a nonlinear diffusion equation for the porosity ϕ, given in Eq. (13). Regardless of the choice of constitutive relations, two boundary conditions and one initial condition are required. At the top of the sample, as we mentioned before, the ice and airflow must go to zero; therefore, the first boundary condition is

$$\begin{array}{}\text{(17)}& w_{\text{i}}=\mathrm{0}\text{at}z=\mathrm{0}.\end{array}$$

In terms of the porosity, this boundary condition translates into

$$\begin{array}{}\text{(18)}& {\displaystyle \frac{\partial \mathit{\varphi }}{\partial z}}=\mathrm{0}\text{at}z=\mathrm{0},\end{array}$$

which is a Neumann boundary condition.

At the bottom of the sample z=h, we apply a constant displacement rate, W; i.e.,

$$\begin{array}{}\text{(19)}& w_{\text{i}}=W\text{at}z=h,\end{array}$$

or equivalently

$$\begin{array}{}\text{(20)}& \left[-N^{\prime }\left(\mathit{\varphi }\right)\right]{\displaystyle \frac{k\left(\mathit{\varphi }\right)}{\mathit{\mu }}}{\displaystyle \frac{\partial \mathit{\varphi }}{\partial z}}=W\text{at}z=h,\end{array}$$

which is also a Neumann boundary condition. Importantly, the constitutive equations play a role at the surface and not at depth. The air pressure at the bottom of the sample is approximately atmospheric, i.e., p=0, due to the ease at which air can escape the sample. The varying load required to compress the sample at a constant rate is

$$\begin{array}{}\text{(21)}& \mathrm{\Sigma }=N\left(\mathit{\varphi }\right)\text{at}z=h,\end{array}$$

which we will compare to the Wang and Baker (2013) experimental data in the next section.

The initial condition is that the snow sample starts with a uniform porosity throughout; i.e.,

$$\begin{array}{}\text{(22)}& \mathit{\varphi }={\mathit{\varphi }}_{\mathrm{0}}\text{at}t=\mathrm{0}.\end{array}$$

Additionally, the height of the snow sample evolves as a free boundary during compaction according to

$$\begin{array}{}\text{(23)}& {\displaystyle \frac{\text{d}h}{\text{d}t}}=-W,\end{array}$$

with the initial condition

$$\begin{array}{}\text{(24)}& h=h_{\mathrm{0}}\text{at}t=\mathrm{0}.\end{array}$$

Figure 3 Compaction profiles for the theoretical model showing solid fraction change (1−ϕ, where ϕ is the porosity) with depth (cf. Hewitt et al., 2016b). Color map shows nondimensional stress and the black lines show the profiles at nondimensional time $t/t_{\mathrm{0}}=\mathrm{0.19},\mathrm{0.38},\mathrm{0.57},\mathrm{0.76},\mathrm{0.92}$. Panels: (a) slow compaction (γ=1000), where the solid fraction is uniform with depth; and (b) fast compaction (γ=1), where significant compaction occurs near the lower part of the sample.

3.3 Nondimensional equations and numerical solution

In our theory for the snow compaction experiments of Wang and Baker (2013), we solve for the evolution of porosity with depth and time, as determined by Eq. (13) with the constitutive relations (14) and (16), the boundary conditions (18) and (20), and the initial conditions (22) and (24). Taking the initial sample height h₀ to be a representative length scale, the displacement rate W to be a characteristic velocity, the initial height to displacement rate as a timescale $t_{\mathrm{0}}=h_{\mathrm{0}}/W$, and k₀ and N₀ to be scales for the permeability and effective pressure, respectively, we can nondimensionalize the variables as $z\to h_{\mathrm{0}}\widehat{z}$, $t\to t_{\mathrm{0}}\widehat{t}$, $k\to k_{\mathrm{0}}\widehat{k}$, $N\to N_{\mathrm{0}}\widehat{N}$, where the hats represent nondimensionalization, although, for ease of notation, we immediately drop the hats. Inserting this nondimensionalization into the equations, we have the following:
the governing equation

$$\begin{array}{}\text{(25)}& {\displaystyle \frac{\partial \mathit{\varphi }}{\partial t}}=\mathit{\gamma }{\displaystyle \frac{\partial }{\partial z}}\left\{(\mathrm{1}-\mathit{\varphi })\left[-N^{\prime }\left(\mathit{\varphi }\right)\right]k\left(\mathit{\varphi }\right){\displaystyle \frac{\partial \mathit{\varphi }}{\partial z}}\right\},\end{array}$$

constitutive equations

$$\begin{array}{}\text{(26)}& N={\displaystyle \frac{(\mathrm{1}-\mathit{\varphi }{)}^n}{{\mathit{\varphi }}^m}}\text{and}k\left(\mathit{\varphi }\right)={\displaystyle \frac{{\mathit{\varphi }}^a}{(\mathrm{1}-\mathit{\varphi }{)}^b}},\end{array}$$

boundary conditions

$$\begin{array}{}\text{(27)}& \begin{array}{rl}& {\displaystyle \frac{\partial \mathit{\varphi }}{\partial z}}=\mathrm{0}\text{on}z=\mathrm{0};\\ & \mathit{\gamma }\left[-N^{\prime }\left(\mathit{\varphi }\right)\right]k\left(\mathit{\varphi }\right){\displaystyle \frac{\partial \mathit{\varphi }}{\partial z}}=\mathrm{1}\text{and}{\displaystyle \frac{\text{d}h}{\text{d}t}}=-\mathrm{1}\text{on}z=h,\end{array}\end{array}$$

and initial conditions

$$\begin{array}{}\text{(28)}& \mathit{\varphi }={\mathit{\varphi }}_{\mathrm{0}}\text{and}h=\mathrm{1}\text{at}t=\mathrm{0},\end{array}$$

where

$$\begin{array}{}\text{(29)}& \mathit{\gamma }={\displaystyle \frac{k_{\mathrm{0}}{N}_{\mathrm{0}}}{\mathit{\mu }{h}_{\mathrm{0}}W}}\end{array}$$

is the ratio of pressure gradient N₀∕h₀ to viscous resistance μW∕k₀. Alternatively, this group can be thought of as the ratio of effective pressure N₀ to viscous pressure μW∕h₀ multiplied by the Darcy number $Da=k_{\mathrm{0}}/h_{\mathrm{0}}^{\mathrm{2}}$ (Bear, 1972).

We solve Eqs. (25)–(28) numerically using a finite volume discretization and the method of lines implemented in Python (LeVeque, 2002). To facilitate the numerical computations, we map the compaction domain into a stationary domain by using the change of variables $\mathit{\xi }=z/h$, which introduces a fictitious advection term into Eq. (25) (Hewitt et al., 2016b). Solutions for two different values of γ and the standard exponents (i.e., a=3, b=2, n=3, m=2) are shown in Fig. 3. The two panels show the two primary regimes: for large values of γ (left panel, γ=1000; Fig. 3a), the porosity ϕ (or equivalently solid fraction, 1−ϕ) is constant with depth. For smaller values of γ (right panel, γ=1; Fig. 3b), compression begins at the bottom of the sample and the top of the sample remains at the starting porosity. As time goes on the sample is compressed and more of the sample compacts.

Figure 4 Comparison of theory with French-press experiments from Wang and Baker (2013): four snow samples were compacted in a Skyscan micro CT scanner at a constant displacement rate (see Sect. 2 and Fig. 1). Using the theory outlined in Sect. 3, the starting density of each sample, and reasonable parameters, we show agreement between theory and experimental data. For the theory, the constitutive exponents are the typical values (a=3, b=2, m=2, n=2), the stress was scaled (i.e., prefactor N₀=30 kPa, friction σ_f=3 kPa), and the values of γ for each curve are shown on the right.

4 Results and discussion

Here we compare the predictions of the theory outlined in the previous section with the experimental data of Wang and Baker (2013), which we described in the Sect. 2. A key output of these experiments was the applied load required to compact the snow as a function of displacement, during the constant-displacement-rate experiments. These data are shown in Fig. 4. The lowest density snow SLT-1 withstands the least stress as a function of displacement. As the density of the snow samples increases, so does the required stress for a given displacement. A similar observation can be made for the theory: the applied load is given by Eq. (21) and can be written dimensionally as

$$\begin{array}{}\text{(30)}& \mathrm{\Sigma }=N_{\mathrm{0}}{\displaystyle \frac{{\left(\mathrm{1}-\mathit{\varphi }\right)}^n}{{\mathit{\varphi }}^m}}\text{at}z=h,\end{array}$$

which states that the applied load Σ increases as the porosity ϕ decreases (see Fig. 2). In other words, for a given displacement, the stress in the theory will also increase as the initial density increases, which is also true for the experimental data, meaning that we have chosen a sensible constitutive law. The theory for the load as a function of the displacement makes sense qualitatively, but we can compare the theory to the experiments in the following way: we run the full model, Eqs. (25)–(28), and evaluate the stress in Eq. (30) using the value of the porosity ϕ at the bottom of the sample z=h and plot it against the displacement $\mathit{\delta }=h_{\mathrm{0}}-h$. The results for reasonable parameters are shown as black lines in Fig. 4.

The model, therefore, contains three inputs and seven parameters. The most important and most uncertain parameter is γ. Due to the fact that the permeability prefactor k₀ and the plastic effective pressure prefactor N₀ are unconstrained, we choose representative values for each sample. The permeability prefactor k₀ only appears within γ, whereas the plasticity prefactor N₀ appears in γ and is required to scale the theory output to compare with experiments. The next four parameters are the exponents a, b, m, and n, which have expected values, but again have not been measured for these samples. Therefore, we only examine small departures from typical values. Additionally, as we will see, there is a small amount of friction in the compression device, and so a small constant stress σ_f is added to the stress output from the theory; i.e., $\mathrm{\Sigma }=N\left(\mathit{\varphi }\right)+{\mathit{\sigma }}_f$ at z=h. This approach was also used by Hewitt et al. (2016b) in their experiments. The inputs are the initial porosity of the sample, which is set by the initial density measured by Wang and Baker (2013). Another input to the problem is the initial height of the sample, h₀. The total experiment run time t_f is also an input and set by the final displacement divided by the displacement rate; i.e., $t_f=(h_{\mathrm{0}}-h_f)/W$, which is about 24 min.

The values of the fixed parameters used for the theory lines in Fig. 4 are the typical exponents (i.e., a=3, b=2, m=2, n=2) and the load information σ_f=3 kPa and N₀=30 kPa. The value of γ for each curve is shown on the figure and is γ=0.46 (SLT-1), γ=0.25 (SLT-2), γ=0.29 (SLT-3), and γ=0.18 (SLT-4). In selecting these values, we successfully (1) found fixed load parameters that worked for all of the experiments and (2) kept the constitutive exponents as the typical values. Thus, γ is the only parameter that varies between theory predictions for each experiment, which is to be expected due to potential variation of k₀ with grain size. It is worth emphasizing that this is not a systemic nonlinear best-fit analysis for seven parameters; rather, the agreement between theory and experiment demonstrates a physically motivated model with a single underconstrained parameter.

In general, there is excellent agreement between the theoretical predictions and the experimental data. For the bottom two samples, SLT-3 and SLT-4, the theory captures all of the major data trends and all falls well within any experimental scatter. For the top sample SLT-1, the theory does a reasonable job in capturing the data trend, yet for these parameters it does not follow the data points exactly. A better fit can be obtained by changing the parameters, indicating that this sample may be a different compaction regime than SLT-3 or SLT-4. However, the fit is reasonable enough that this sample is likely just on the edge of a new regime, if at all. In contrast to the other samples, the theory does not adequately capture the data trend of the middle sample, SLT-2. This sample is interesting because it has almost the same density and specific surface area as SLT-3, yet it responded very differently to compression loading. Due to the small size of the samples, i.e., 15.7 mm in diameter and 18 mm tall, the likelihood of defects or inhomogeneities dominating the results is quite high. Thus, the snow bonds in SLT-2 from prior sintering could have been particularly stubborn, requiring more load for a given displacement.

Another possibility for why the two snow samples SLT-1 and SLT-2 did not agree as well with the theory as SLT-3 and SLT-4 is that pressure sintering during the experiment allowed for greater bonding of snow crystals. Adding pressure is an efficient method of accelerating the rate of sintering and can lead to sinter rates that are orders of magnitude faster than by ambient surface energy differences alone (Rahaman, 2007). For this reason, Wang and Baker (2013) attribute the increase in required load to accelerating the sintering and coarsening processes occurring within the snow samples. Willibald et al. (2019) also analyze sintering during compaction experiments and find that the sintering rate is enhanced. The plastic compaction theory we present in this paper does not include pressure processes and, therefore, would not be able to describe the stress required to break sintered bonds, although this is a promising direction for future research.

The plastic compaction theory presented here can be related back to the general firn compaction model given in Eq. (1). Rearranging Eq. (6) gives

$$\begin{array}{}\text{(31)}& {\displaystyle \frac{\partial \mathit{\varphi }}{\partial t}}+w_{\text{i}}{\displaystyle \frac{\partial \mathit{\varphi }}{\partial z}}=(\mathrm{1}-\mathit{\varphi }){\displaystyle \frac{\partial w_{\text{i}}}{\partial z}},\end{array}$$

which is in the form of Eq. (1). Inserting Darcy's law (9) with mass conservation (5) and (6), the compaction function 𝒞 is given by

$$\begin{array}{}\text{(32)}& \begin{array}{rl}\mathcal{C}& =-\left({\mathit{\rho }}_{\text{a}}-{\mathit{\rho }}_{\text{i}}\right)\left(\mathrm{1}-\mathit{\varphi }\right){\displaystyle \frac{\partial w_{\text{i}}}{\partial z}}\\ & =\left({\mathit{\rho }}_{\text{a}}-{\mathit{\rho }}_{\text{i}}\right)\left(\mathrm{1}-\mathit{\varphi }\right){\displaystyle \frac{\partial }{\partial z}}\left[N^{\prime }\left(\mathit{\varphi }\right){\displaystyle \frac{k\left(\mathit{\varphi }\right)}{\mathit{\mu }}}{\displaystyle \frac{\partial \mathit{\varphi }}{\partial z}}\right],\end{array}\end{array}$$

which shows the connection between compaction and air evacuating pore space as well as the role of the constitutive relations N(ϕ) and k(ϕ). In this way, measuring the parameters of these constitutive relations in the laboratory allows for predictions of compaction using 𝒞 from Eq. (32).

For generalized viscous compaction, McKenzie (1984) starts from Eq. (31) and connects the divergence $\partial w_{\text{i}}/\partial z$ to an effective pressure P_e through a compaction viscosity η_c; i.e.,

$$\begin{array}{}\text{(33)}& {\displaystyle \frac{\partial w_{\text{i}}}{\partial z}}={\displaystyle \frac{P_{\mathrm{e}}}{{\mathit{\eta }}_{\mathrm{c}}}},\end{array}$$

which is the compaction law used in studies of temperate ice (Schoof and Hewitt, 2016; Hewitt and Schoof, 2017; Meyer et al., 2018) and is equivalent to the viscous closure of a Röthlisberger channel (Nye, 1953; Fowler, 1984; Meyer et al., 2016, 2017). The McKenzie (1984) compaction law, Eq. (33), implies that

$$\begin{array}{}\text{(34)}& \mathcal{C}=-\left({\mathit{\rho }}_{\text{a}}-{\mathit{\rho }}_{\text{i}}\right)\left(\mathrm{1}-\mathit{\varphi }\right){\displaystyle \frac{P_{\mathrm{e}}}{{\mathit{\eta }}_{\mathrm{c}}}}.\end{array}$$

If we take the effective pressure P_e to be the overburden pressure, i.e., P_e=ρgz, thereby setting the air pressure to zero, we find the pressure-compaction model described by Cuffey and Paterson (2010) following Arthern et al. (2010). Although this compaction model is attractive because it connects viscous processes to compaction, by setting the air pressure to zero (or even a constant), this compaction model fails to capture the evacuation of air from the pore space, which we have shown to be a very important process in the mechanical compaction of laboratory snow samples. In future work, we will adapt the model we present here to include viscous compaction following McKenzie (1984). Instead of neglecting the air pressure, as in the Arthern et al. (2010) work, we will explicitly solve Darcy's law and determine the resulting compaction.

5 Conclusions

In this paper, we articulated a mathematical model to describe the snow compaction experiments of Wang and Baker (2013). This model consists of mass and momentum conservation as well as porosity-dependent permeability and plastic effective pressure constitutive equations. The outputs of the model are the snow density as a function of time and space as well as the stress as a function of displacement. Comparing the model outputs to the experimental data of Wang and Baker (2013) shows excellent agreement, especially for the low-density sintered snow. As the density increased, small discrepancies between model and theory emerged, potentially due to the necessity of creep or sample defects. Nevertheless, the excellent agreement between theory and experiments suggests that measuring compaction in the lab is a promising direction forward for understanding snow compaction and that the plastic effective pressure as a function of porosity may be a key constitutive relation to quantify.