Introduction

The densification of snow has been studied by a number of investigators, including the pioneering works of Reference BaderBader (1960, Reference Bader and Kingery1963), Reference BensonBenson (1962), and Reference Anderson, Benson and KingeryAnderson and Benson (1963). These authors calculated the density of dry snow for different depths. A study of the age of dry snow as a function of depth was done recently by Reference Herron and LangwayHerron and Langway (1980), with good results, by treating the snow as two layers.

In this paper we try to answer the questions: (1) Can an age–depth relation for dry snow as one continuous layer be derived without use of an empirical formula for the density-time relation? (2) What kind of stress–strain relation is required to arrive at the steady-state density profile described by Reference LingLing (1985)?

Let us first consider snow at constant accumulation and constant temperature. Then, the process of snow densification is invariant with respect to the snow surface (Reference BaderBader, 1960). If we want to know what the density of the snow-surface layer would be 10 years from now, all we have to do is to dig into the snow-pack and find the layer that is 10 years old and measure its density.

The Stress–Strain Relation

Consider the steady-state profile equation for dry snow from Reference LingLing (1985), where z is measured positive downward from the surface,

This is based on a non-linear differential equation relation between the change of pressure p and the change of density ρ:

where ρ is the density at depth z, ρ ₀ is the surface density of snow, ρ _m is the maximum attainable density, and L is a characteristic length scale, set equal to 38 m according to Reference LingLing (1985). Equation (1a) is obtained by (1) generalization of a simple physical model dρ = c(ρ _m − ρ)dρ to the form d(ρ ⁿ) = c(ρ _m − ρ)dp, and (2) by setting n = 2 for mathematical convenience, as in Reference LingLing (1985). Differentiating Equation (1) with respect to z yields

For any snow layer of infinitesimal thickness dz, the mass per unit area is

where A is the accumulation rate in m/year of water equivalent, ρ _w is the density of water, and dt is the infinitesimal difference between the time a particle at the bottom of the layer was deposited and the time a particle at the top was deposited. Using Equations (1) and (3), and integrating with the assumption of constant accumulation rate A, one gets

Equation (4) gives snow age as a function of snow depth, for constant temperature. Combining Equations (1), (2), and (3) to eliminate z yields

Integrating Equation (5), from ρ = ρ ₀ at time = 0 yields

Now the pressure at a layer, taking A to be constant in time, is the total weight per unit area above this layer, which can be obtained by multiplying Equation (3) by g and then integrating from t = 0 to t = t:

where t = 0 is the time when the particle was at the surface. Substituting Equation (7) into Equation (5), to relate density to pressure, yields the stress–strain relation

This indicates that the deformation rate of dry snow is proportional to the pressure or overload as well as to the factor (ρ _m − ρ) and is inversely proportional to the time t When t = 0, the deformation rate will be infinite unless ρ ∞ t also. That a sudden load on a layer of snow causes infinite strain has been shown from results of work done by Reference Colbeck, Colbeck, Shaw and LemieuxColbeck and others (1978), and Reference Costes and KingeryCostes (1963). Therefore, such behavior in the equation is quite desirable and realistic.

Substituting Equation (6) into Equation (8) to eliminate t now gives

This is the stress–strain relation for dry snow which we will call a “non-linear stress–strain relation” since it is derived from d(ρ ²) = c(ρ _m − ρ) dp. However, the usual way of writing the stress–strain relation for dry snow, as used by Reference BaderBader (1960), Reference MellorMellor (1964), and Reference AndersonAnderson (1976), among others, is dρ/ρdt = p/η where η is a “compactive viscosity factor”. Reference BaderBader (1960) pointed out that the assumption of Newtonian viscosity as written above can be used when p < 1000 g/cm², which occurs at depths of between 15 and 20 m. Here we will use it for depth far beyond this range to see how it works. Equation (9) thus gives

Re-arranging Equation (9) to integrate from ρ = ρ ₀ at time t = 0, gives

If the density ratio r(ρ) = ρ/ρ _m is introduced, from which dρ = ρ _mdr, the first integral of Equation (10) may be written as

in which r ₀ indicates r(ρ ₀). If the change of variable θ(r) = (1 − r) − ln(1 − r) is made, from which dr = (1 − r) dθ/r, it becomes

which can be directly integrated to give

The inverse function r(F) must be obtained if ρ is to be determined from a known value of the second integral of Equation (10), which is equal to F. Thus, given a value of F, first the inverse of Equation (11) is used to get r, and then ρ is recovered from r by using the relation r = ρ/ρ _m The dependence of r on F is shown in Figure 1a for selected values of the parameter r ₀.

Fig. 1.

a. Inverse of density integral, given by Equation (10), for indicated values of the parameter r₀. The curves are tangent to the F = 0 axis at r = r₀ and they asymptotically approach r = 1 as F → ∞. b. Error curves, for indicated values of the parameter r₀, arising from use of Equation (11) to approximate the inverse of the density integral. Each curve, whose coefficients are given in Table I, minimizes the maximum error |r* − r| over the interval r₀ ≤ r ≤ 1, at each end of which the error is identically zero.

Table I.

The Optimizing Coefficients of the f-Inverse Approximation, Equation (11)

Approximating the Inverse of the Density Integral

Because an inverse of Equation (11) has not been found, it is not possible to get r directly from any particular value of F. Instead, a simple empirical approximation is devised to get an estimate r ^* of it, given a value of F. It is constructed simply by matching a mathematical function closely to points on the inverse function r(F).

While constructing the approximation, points r(F) are obtained numerically. As Equation (11) readily gives F for any value of r, it is possible to search for that value of r corresponding to a particular desired F. In practice, it would be tedious to undertake such an interation any time one wanted to get r from F, so the approximating function is constructed as a reasonably accurate and very convenient way of estimating it.

The form chosen here for r ^*(F) has exactly the correct behavior at the ends of the interval: r ^* = r ₀ when F = 0, and r ^* = 1 when F → ∞. The normalization on ρ _m is extended here by considering the ratio F/ρ ² _m instead of F itself. The devised function is

The empirically determined values of the two coefficients a and b are chosen to minimize the maximum value of |r ^* − r| occurring over the interval r ₀ ≤ r ≤ 1. Just as r ₀ is a parameter of the function F(r), it is also a parameter of its approximation. Thus, the minimizing values of a and b depend on the value of r ₀ (Table I). The error r ^* − r is shown as a function of r in Figure lb for selected values of the parameter r ₀. Other functional forms could be used to approximate the inverse function more accurately, but at the cost of greater algebraic and computational complexity.

The accuracy of the approximation would be improved if the requirements of exact behavior at the end points were relaxed. For the case r ₀ = 0.1, relaxing the r ^* = 1 condition at F = ∞, by replacing the factor 1 − ρ ₀ by a third free coefficient, would reduce the maximum error |r* − r| by 8%; relaxing the r ^* = r ₀ condition at F = 0 as well, by replacing the r ₀ term by a fourth free coefficient, would reduce it by an additional 18%. Also for the r ₀ = 0.1 case, shortening to r ₀ ≤ r ≤ r _max the interval over which the maximum error |r ^* − r| is to be minimized, by using different values for the two coefficients a and b, would reduce the maximum error by 16% for r _max = 0.7, by 3% for 0.8, and by none for 0.85 or greater. Because the error is already of the order of density-measuring error, and because achieving those small error reductions would impair the convenience of using the approximation, the values of the coefficients given in Table I are used here with the simple form given by Equation (12); the end-point conditions are met, and the maximum error is minimized over the full interval r ₀ ≤ r ≤ 1.

The Temperature Correction

Equation (10) is for dry snow at constant temperature; but the effects of temperature variation on the densification of snow can be quite substantial, and therefore should be included. The snow temperature varies annually in the upper 10 m so that values significantly higher than the mean value, T _m, are encountered each year. Reference Bader and KingeryBader (1963) corrected for this by calculating an equivalent temperature, T _e, which is higher than T _m because the higher temperatures in the annual variations are most effective in modifying the rates of densification in the top 10 m. We use the temperature factor β = β ₁exp(−E/RT) as in Reference GlenGlen (1955), Reference Bader and KingeryBader (1963), Reference AndersonAnderson (1976), and Reference Herron and LangwayHerron and Langway (1980). One may write (1/ρ)(dρ/dt) ∞ β ₁expt(−E/RT) but (1/ρ)(dρ/dt) ∞ p as in Equation (9). Thus (1/ρ)(dρ/dt) ∞ p β ₁exp(−E/RT). To include the temperature effect, we need only multiply p by β ₁exp(−E/RT) on the right-hand side of Equation (10). When T = T _m, β = 1, so β ₁ = exp(E/RT _m) and β = exp(E/RT _m − E/RT). Here E is the activation energy, which is taken to be 1.33 × 10⁵J/mole according to Reference GlenGlen (1955). The gas constant R is taken to be 8.3 J/kmole. T is the temperature of the snow, and T _m is the temperature at the depth of 10 m. At this depth the annual variation in temperature is less than 0.5 deg.

We assume a temperature distribution as follows: based on the work of Reference Weller and SchwerdtfegerWeller and Schwerdtfeger (1968), and Reference BensonBenson (1962), where ΔT is the annual temperature amplitude (that is, half the total annual range) at the snow surface of a specific station, assumed to be 15 deg in this work, and τ is the period, taken here to be 1 year, and α is the thermal diffusivity which is equal to k/ρc where k is the thermal conductivity of snow and c is the specific heat. The thermal conductivity of snow varies from 0.04 to 1.22 J m⁻¹ s⁻¹ for snow up to 2.5 m in depth, as shown by Reference LangeLange (1985) in his work on measurements in the Antarctic. Here, we use a representative Value of k = 0.837 J m⁻¹ s⁻¹ K⁻¹ for the upper 10 m of snow. Now, ρ = 0.4 Mg m⁻³ and c = 1.967 × 10⁻³ J kg⁻¹ K⁻¹ (from Reference MellorMellor, 1977), so α = 1.064 × 10⁻⁶ m² s⁻¹. Thus, Equation (10), along with Equation (11), with the inclusion of the temperature factor, is

where T is as given above. The t coordinate is chosen so that t = 0 in the middle of the year when the surface layer was deposited, and this is done uniformly for any specific layer during the calculation. Depth (z) is always available as a measured quantity. However, since we need to integrate Equation (13), we need depth in terms of time. This can be approximated by using Equation (4). After expanding e^−(z/L) into a Taylor series and keeping the first three terms:

Therefore,

Equation (12) is used to get ρ from the right-hand side of Equation (13) and the corresponding depth, from Equation (3), is

where the density distribution as a function of t is obtained from Equation (13).

Results

Equations (13) and (14) have been used to calculate the age of snow layers at five stations — Crête, Site 2, and Milcent of Greenland; and Byrd Station and Little America V, Antarctica. Table II shows the input data used for the five stations. The relations between observed age and depth for the five stations are from Reference Herron and LangwayHerron and Langway (1980), and Reference GowGow (1968). Figure 2 shows comparisons of calculated age at given depth with observed data, with and without the temperature correction. Figure 3 shows comparisons of calculated density–depth curves with observation for Byrd Station and Little America V (the two stations for which tabulated density–depth data are available to the authors), with and without the temperature correction. There does not seem to be much difference for the age versus depth calculation where we include the temperature correction; however, there is some improvement for the calculation of the density versus depth when we include the temperature correction, as shown in Figure 3, even though for Byrd Station the agreement is still not very satisfactory in the upper 20 m (Fig. 3a).

Table II.Parameters Used for Five Snow Stations

Fig. 2.

Snow age versus depth for selected sites. Full curves show calculation with temperature variation, and dashed curves show calculation with constant temperature. Actual values (x) are from Reference Herron and LangwayHerron and Langway (1980).

Fig. 3.

Snow density versus depth for selected sites. Full curves show calculation with temperature variation, and dashed curves show calculation with constant temperature. Actual values (x) are from Reference Bader and KingeryBader (1963) and Reference GowGow (1968).

The lack of agreement in the upper 20 m results because we have ignored the discontinuity in physical properties which occurs at a porosity of about 40%, i.e. a density of 0.55 Mg m⁻³ (Reference BensonBenson, 1962; Reference Anderson, Benson and KingeryAnderson and Benson, 1963). The significance of the change in predominant densification mechanisms at this “critical density” in Antarctic snow was emphasized by Reference Robertson and BentleyRobertson and Bentley (1975) in their study of seismic velocity gradients. Equation (12) was also used in the above calculations; the approximations are so close to the exact calculations that the differences are negligible (see Table I for the maximum error). Thus, the approximation in Equation (12) can be used with confidence in actual calculations.

That the results of the calculation are good leads us to believe that the stress-strain relation, Equations (9) and (10), is usable for Greenland and Antarctic snow.

The accuracy of the results is comparable to that of Reference Herron and LangwayHerron and Langway (1980). Our model is a one-layer continuous formulation, while the model by Herron and Langway is discontinuous at the density of 0.55 Mg m⁻³.

Conclusion

We have developed a stress–strain relation for dry snow in Greenland and Antarctica. It is used to calculate the age–depth relation of dry snow, and the results are quite encouraging when given the accumulation rate, surface-snow density, and deep-core temperature. The average error with (Fig. 3) respect to the observed age is less than 3%, while the maximum error is about 7.5%.

With temperature correction for the top 10 m of snow, we are able to improve on density–depth calculations.

We believe this method will be valid only in the case of dry snow in the cold regions. However, it could be extended for dry snow in other regions by changing parameters such as the characteristic length L.

When one wants to study wet snow or snow under strong thermodynamic influences such as melting and freezing, perhaps the method has to be modified or a totally different approach might be needed.