The densification of snow has been studied by a number of investigators, including the pioneering works of Bader (1960, 1963), Benson (1962), and Anderson 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 Herron 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 Ling (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 (Bader, 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 Ling (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, ρ
0 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 Ling (1985). Equation (1a) is obtained by (1) generalization of a simple physical model dρ = c(ρ
m − ρ)dρ to the form d(ρ
n) = c(ρ
m − ρ)dp, and (2) by setting n = 2 for mathematical convenience, as in Ling (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 ρ = ρ
0 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 Colbeck and others (1978), and Costes (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(ρ
2) = c(ρ
m − ρ) dp. However, the usual way of writing the stress–strain relation for dry snow, as used by Bader (1960), Mellor (1964), and Anderson (1976), among others, is dρ/ρdt = p/η where η is a “compactive viscosity factor”. Bader (1960) pointed out that the assumption of Newtonian viscosity as written above can be used when p < 1000 g/cm2, 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 ρ = ρ
0 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
0 indicates r(ρ
0). 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 r0. The curves are tangent to the F = 0 axis at r = r0 and they asymptotically approach r = 1 as F → ∞. b. Error curves, for indicated values of the parameter r0, 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 r0 ≤ 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
0 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
0 ≤ r ≤ 1. Just as r
0 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
0 (Table I). The error r
* − r is shown as a function of r in Figure lb for selected values of the parameter r
0. 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 = 0.1, relaxing the r
* = 1 condition at F = ∞, by replacing the factor 1 − ρ
0 by a third free coefficient, would reduce the maximum error |r* − r| by 8%; relaxing the r
* = r
0 condition at F = 0 as well, by replacing the r
0 term by a fourth free coefficient, would reduce it by an additional 18%. Also for the r
0 = 0.1 case, shortening to r
0 ≤ 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
0 ≤ 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. Bader (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 β = β
1exp(−E/RT) as in Glen (1955), Bader (1963), Anderson (1976), and Herron and Langway (1980). One may write (1/ρ)(dρ/dt) ∞ β
1expt(−E/RT) but (1/ρ)(dρ/dt) ∞ p as in Equation (9). Thus (1/ρ)(dρ/dt) ∞ p
1exp(−E/RT). To include the temperature effect, we need only multiply p by β
1exp(−E/RT) on the right-hand side of Equation (10). When T = T
m, β = 1, so β
1 = exp(E/RT
m) and β = exp(E/RT
m − E/RT). Here E is the activation energy, which is taken to be 1.33 × 105J/mole according to Glen (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 Weller and Schwerdtfeger (1968), and Benson (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−1 s−1 for snow up to 2.5 m in depth, as shown by Lange (1985) in his work on measurements in the Antarctic. Here, we use a representative Value of k = 0.837 J m−1 s−1 K−1 for the upper 10 m of snow. Now, ρ = 0.4 Mg m−3 and c = 1.967 × 10−3 J kg−1 K−1 (from Mellor, 1977), so α = 1.064 × 10−6 m2 s−1. 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:
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).
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 Herron and Langway (1980), and Gow (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 Herron 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 Bader (1963) and Gow (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−3 (Benson, 1962; Anderson and Benson, 1963). The significance of the change in predominant densification mechanisms at this “critical density” in Antarctic snow was emphasized by Robertson 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 Herron 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−3.
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.
We want to thank A. Thorndike for his suggestions. This work would not have been done without his inspiring discussions. We also want to thank S. Hodge, E.G. Josberger, and W.J. Campbell for their suggestions. The help from R. Armstrong is also appreciated.
Anderson, D.L., and
1963. The densification and diagenesis of snow. (In Kingery, W.D., ed. Ice and Snow: Properties. Processes, and Applications; proceedings of a conference held at the Massachusetts Institute of Technology, February 12–16, 1962. Cambridge, MA, M.I.T. Press, p. 391–411.)
1976. A point energy and mass balance model of a snow cover. Silver Spring, MD, National Weather Service, Office of Hydrology. (PB-254 653.)
1960. Theory of densification of dry snow on high polar glaciers. U.S. Army Snow. Ice and Permafrost Research Establishment. Research Report 69.
1963. Theory of densification of dry snow on high polar glaciers II. (In Kingery, W.D., ed. Ice and Snow; Properties, Processes, and Applications; proceedings of a conference held at the Massachusetts Institute of Technology. February 12–16, 1962. Cambridge, MA, M.I.T. Press p. 351–76.)
1962. Stratigraphic studies in the snow and firn of the Greenland ice sheet. U.S. Army Snow, Ice and Permafrost Research Establishment. Research Report 70.
Colbeck, S.C., and others. 1978. The compression of wet snow, Colbeck, S.C., Shaw, K.A., Lemieux, G.
1963. Confined compression tests in dry snow. (In Kingery, W.D., ed. Ice and Snow; Properties, Processes, and Applications; proceedings of a conference held at the Massachusetts Institute of Technology. February 12–16, 1962. Cambridge, MA, M.I.T. Press p. 594–612..)
1955. The creep of polycrystalline ice. Proceedings of the Royal Society, Ser. A, Vol. 228, No. 1175, p. 519–38.
1968. Deep core studies of the accumulation and densification of snow at Byrd Station and Little America V, Antarctica. CRREL Research Report
Herron, M.M., and Langway, C.C. jr. 1980. Firn densification: an empirical model. Journal of Glaciology, Vol. 25, No. 93, p. 373–85.
1985. Measurements of thermal parameters in Antarctic snow and firn. Annals of Glaciology, Vol. 6 p. 100–04.
1985. A note on the density distribution of dry snow. Journal of Glaciology, Vol. 37 No. 108, p. 194–95.
1977. Engineering properties of snow. Journal of Glaciology, Vol. 19 No. 81, p. 15–56.
Robertson, J.D., and Bentley, C.R.
1975. Investigations of polar snow using seismic velocity gradients. Journal of Glaciology, Vol. 14 No. 70, p. 39–48.
Weller, G., and Schwerdtfeger, P.
1968. Thermal properties and heat transfer processes of the snow of the central Antarctic plateau. [Union Géodésique el Géophysique Internationale. Association Internationale d’Hydrologie Scientifique.] [International Council of Scientific Unions. Scientific Committee on Antarctic Research. International Association of Scientific Hydrology. Commission of Snow and Ice.] International Symposium on Antarctic Glaciological Exploration (ISAGE), Hanover, New Hampshire, U.S.A., 3–7 September 1968, p. 284–98 [(Publication No. 86 [de l’Association Internationale d’Hydrologie Scientifique].)].