Introduction

The transformation of snow into glacial ice is one of the most fundamental processes in glaciology. A number of theories regarding the dry densification of snow and firn have been proposed Reference Robin(Robin, 1958; Reference BensonBenson, 1959, Reference Benson1962; Reference BaderBader, 1960, Reference Bader1962; Reference Costes and KingeryCostes, 1963; Reference Anderson, Benson and KingeryAnderson and Benson, 1963; Reference Kojima and MellorKojima, 1964; Reference HobbsHobbs, 1968; Reference GowGow, 1975; Reference JohnsenJohnsen, 1977) and yet the relative importance of possible component processes such as settling, sublimation, recrystallization, volume diffusion, surface diffusion, and other deformation processes has not been firmly established. Consequently, it has not been possible to develop a theoretically based model which accurately predicts the depth–density profile from a given site without first obtaining seismic compressional-wave velocity data Reference Kohnen(Kohnen, 1972) or laboratory measurements on deep firn samples from that location.

Three stages of densification have been identified. The densification rate in the first stage, from initial snow density to the “critical density” of about 0.55 Mg m⁻³ Reference Benson(Benson, 1962; Reference Anderson, Benson and KingeryAnderson and Benson, 1963), is usually the most rapid. The dominant mechanism is considered to be grain settling and packing. Densities increase more slowly with depth in the second stage until the intercommunicating air passages become closed off to form individual bubbles at densities of 0.82–0.84 Mg m⁻³. The mechanisms involved are not well understood. In the third stage, below the pore close-off zone, air can no longer be excluded and further densification takes place by compression of the bubbles Reference Langway(Langway, 1958; Reference BaderBader, 1965). It is generally recognized that the form of depth–density profiles depends on the snow accumulation rate and temperature. Higher temperatures and lesser accumulation rates result in the greatest change of density with depth.

The purpose of the present work is to derive an empirical, but widely applicable, model of the first and second stages of polar snow densification with the aid of depth-density data from recently recovered Greenland and Antarctic ice cores. Such a model would be of use to glaciologists and ice drilling engineers if the depth–density, depth–load, and depth–age profiles and the approximate pore close-off depth could be predicted from the commonly measured variable of mean annual temperature and mean annual accumulation rate. It is hoped that the model and the density data may also aid others in developing a better physical understanding of the densification process.

The model

Following Sorge’s Law Reference Bader(Bader, 1954), the relation between snow density ρ and the depth h below the snow surface is assumed invariant with time under conditions of constant snow accumulation and temperature. Thinning of annual layers by plastic flow Reference Nye(Nye, 1963), a relatively minor effect in the upper layers of most inland ice sheets, is ignored.

Depth–density information was obtained from 17 sites in Greenland and Antarctica where mean annual temperatures range from −57°C to −15°C and annual accumulation rates range from 0.022–0.5 m water year⁻¹ (Table 1). The drilling techniques and reliability of the data are found in some of the original references. Temperatures used in this model are those at the 10 m depth and measurement accuracy is generally ±0.1 deg or better. All densities were determined from measurements of the mass and dimensions of sections of core. Because of firn permeability, the highly accurate method of hydrostatic weighing Reference Butkovich(Butkovich, 1953) is impractical. Measurement accuracy is frequently ±0.001 Mg m⁻³ although individual data points may exhibit more scatter due to seasonal variations in density or the presence of melt features. There is often a large uncertainty in the annual accumulation values due to imprecision in the identification of annual layers, short-range geographical and temporal variations, and the relatively short time period for which data are available.

The model is based on Reference RobinRobin’s (1958) suggestions that in the densification of firn the proportional change in air space is linearly related to the change in stress due to the weight of overlying snow. Reference SchyttSchytt (1958) expressed this as:

where ρ _i is the density of ice (0.917 Mg m⁻³), which implies a linear relationship between ln [ρ/(ρ _i − ρ)] and depth. Schytt found such a relationship in the Maudhehn core data.

To explore the widespread applicability of Equation (1) and the nature of the “constant” therein, depth-density information from the other available sites was plotted as ln [ρ/(ρ _i − ρ)] versus depth (Appendix). Temperature and pressure corrections for the density of pure ice make insignificant changes in ln [ρ/(ρ _i − ρ)] over the density range involved and a value of ρ _i = 0.917 Mg m⁻³ was used throughout. The plots generally consisted of two linear segments, one for ρ < 0.55 Mg m⁻³, and a steeper slope for 0.55 Mg m⁻³ < ρ < 0.8 Mg m⁻³, corresponding to the first and second stages of densification. Although pore close-off occurs at ρ = 0.82−0.84 Mg m⁻³, many of the plots were linear only through ρ = 0.8 Mg m⁻³ with a much greater densification rate at the zone of pore close off. The pronounced influence of ice layers in cores from warmer sites such as Maudheim Reference Schytt(Schytt, 1958), Roi Baudouin Reference Tongiorgi, Tongiorgi, Picciotto, de Breuck, Norling, Giot and Pantanetti(Tongiorgi and others, 1962), and Dye 2 is evidenced by the scatter in the profiles. Only the ρ < 0.55 Mg m⁻³ portion of the Roi Baudouin site was used in the derivation of the model.

The slopes of the linear segments, fitted by a least-squares method or, in some cases, by eye, may be represented as

and

where C and C′ are constants for each site. To solve for the densification rate, dρ/dt, Equations (2a) and (2b) are rearranged with the substitution dh/dt = A/ρ, where A is the accumulation rate in water-equivalent units:

and

Equations (3a) and (3b) are equivalent to Equation (1). However, since the accumulation rate and temperature are fairly constant at each site, the functional dependencies of C and C′ on T and A are unknown. We have assumed that the temperature and accumulation-rate dependencies may be separated and that the rate equations are of the form

and

where k ₀ and k ₁ are Arrhenius-type rate constants dependent only on temperature and a and b are constants dependent on the densification mechanisms. Values of a and b may be determined by comparing slopes for each stage of densification at sites of nearly equivalent temperature and different accumulation rates using

and similarly for b by substituting C′ for C into Equation (5). Values determined in this manner (Table II.) are a = 1.1±0.2 and b = 0.5±0.2. Taking the values of a and b as 1 and 0.5, respectively, this model agrees with Reference RobinRobin’s (1958) description of the densification process for the first stage, but not for the second.

Values of the rate constants were deduced for each site and the general equations of k ₀ and k ₁ were determined from Arrhenius plots of In k against 1/T (Fig. 1).

and

where R is the gas constant (8.314 J K^{− 1} mol⁻¹) and T is the temperature in kelvins. The correlation coefficients between the logarithms of the rate constants and reciprocal temperature are 0.959 for k ₀ and 0.987 for k ₁. S. J. Johnsen has kindly pointed out that Equation (4b) may be rewritten in terms of load σ as:

for which the apparent activation energy is 42.6 kJ/mol. This value is more in line with values for known physical processes in ice. However, the utility of Equation (4b) lies in the case of calculation of density profiles or accumulation rates when one of these and the mean annual temperature are available.

Fig. 1. Arrhenius plots for k_o (lower line) and k₁ (upper line ).

Depth–density and Depth–Age Calculations

Given T, A, and the initial snow density ρ ₀ at a given site, it is possible to calculate the density at depth, the depth of the ρ = 0.55 Mg m⁻³ transition, and the depth-age relationship down to the zone of pore close-off. The value of ρ ₀ used is the zero-depth intersection of the line of ln [ρ/(ρ _i − ρ)] versus depth, which may not be identical with the commonly reported average density over the first one or two meters of snow.

For the initial stage of densification, the density at depth h is:

where

. Note the surprising result that the depth-density relationship in this region is independent of the accumulation rate. The depth at which the “critical density” is reached is, by rearrangement,

and the age in years at this depth is

Intermediate depths and ages may be calculated by substituting appropriate values of ρ into Equations (8) and (9).

For the second stage of densification, the density at depth h is

where

. The age of firn at density ρ is

The mean annual accumulation rate may be estimated from the slope C’ of the second stage of densification and the mean annual temperature:

The effects of independently varying the temperature and accumulation rate on the profile of ln [ρ/(ρ _i−ρ)] predicted by this model are shown in Figure (2). In Figure (2a), the accumulation rate is held constant at 0.3 m water year⁻¹ and the initial density is 0.36 Mg m⁻³. For temperatures from −15°C to −40°C, the depth at which the “critical density” of 0.55 Mg m⁻³ is reached varies from 9.3 m to 15.5 m. At −15°C the density attains a value of 0.8 Mg m⁻³ at a depth of 44 m and an age of 93 years. The corresponding values at −40°C are 115 m and 254 years.

Fig. 2. Predicted depth—ln [ρ/(ρ_i−ρ)] curves for: (a) various temperatures at a constant accumulation rate of 0.30 m water year⁻¹ and initial density of 0.36 Mg m⁻³, and (b) various accumulation rates (in m water year⁻¹) at a constant temperature of −30°C and initial density 0.36 Mg m⁻³.

In Figure (2b) the temperature is held constant at −30°C, the initial density is again 0.36 Mg m⁻³, and the accumulation rate is varied from 0.1 to 0.6 m water year⁻¹. As noted previously, the first-stage depth-density curve is independent of the accumulation rate. The density reaches a value of 0.8 Mg m⁻³ at 49 m for A = 0.1 compared to a depth of 102 m for A = 0.6 m water year⁻¹. However, this process requires only 113 years for A = 0.6 compared to 310 years for A = 0.1 m water year⁻¹.

Comparison with Observation

Depth–density data from inland Greenland (Site 2, South Dome, North Central), inland Antarctica (Vostok, Old Byrd, Roosevelt Island Dome) and the Ross Ice Shelf, Antarctica (Little America V, J-9, C-7-3) are shown in Figure(3). Model curves, shown as solid lines, were generated from Equations (7) and (10) using the initial density, temperature, and accumulation rates given. The agreement between observed and predicted profiles for sites with such a wide range of temperatures and accumulation rates demonstrates the usefulness of the model. The density anomalies at 22 m in the Old Byrd profile and at 28 m in the Little America V profile will be discussed more fully in the next section.

Fig. 3. Depth-density curves for: (a) Greenland, (b) inland Antarctic, and (c) Ross Ice Shelf ice cores. Initial density, accumulation rate, and temperature values shown were used to generate model curves shown as solid lines. Arrows in the Site 2 profile refer to pronounced melt features Reference Langway(Langway, 1970).

As a further measure of the scope of this model, Table III compares ages of five ice cores calculated from Equations (9) and (11) with ages determined by stratigraphic techniques Reference Gow(Gow, 1968; Reference LangwayLangway, 1970) and by analysis of stable oxygen isotope and microparticle concentration variations Reference Hammer, Hammer, Clausen, Dansgaard, Gundestrup, Johnsen and Reeh(Hammer and others, 1978; Reference Reeh, Reeh, Clausen, Dansgaard, Gundestrup, Hammer and JohnsenReeh and others, 1978; personal communication from W. Dansgaard). Observed and predicted ages agree to within five years over time spans as great as 256 years, which is very good considering the uncertainties in measurement accuracy, the normal variability in accumulation rates, and the empirical nature of the model.

Table IV compares observed accumulation rates with rates derived using Equation (12). Due to the abundance of ice layers at Dye 2, Maudheim, and Roi Baudouin, these sites were not included in the derivation of the model, and predicted values are included for completeness. The agreement between observed and predicted values is quite good with an average deviation of only 0.04 m water year⁻¹ and an average relative deviation of 16%. For dry-snow and percolation facies, the accumulation rales predicted by this model will generally be within about 20% of the mean over the time period represented in the ρ < 0.8 Mg m⁻³ section of core.

Density Profiles as Climatic Indicators

Discrepancies between observed and predicted density profiles should be indicative of climatic conditions, since the pertinent variables are temperature, accumulation rate, and the presence of ice layers. This concept is best demonstrated with the ice cores from Site 2, Old Byrd Station, and Little America V, where long-term accumulation-rate records are available, the density values have been averaged over essentially every meter of depth, and the cores were obtained by mechanical, rather than thermal, drilling techniques Reference Gow(Gow, 1968; Reference LangwayLangway, 1970).

The most pronounced melt features in the 0−70 m depth interval of the Site 2 core Reference Langway(Langway, 1970) are indicated by arrows in Figure (3a). The observed and predicted densities are in excellent agreement throughout most of the core, and a clear correspondence may be seen between the presence of ice layers and abnormally high densities. Without the predicted profile for comparison, such features might have been ascribed to scatter in the data.

Anomalous intervals of constant density may be clearly seen at a depth of 22 m in the Old Byrd Station core and at the 28 m depth in the Little America V core. Above these depths, observed and predicted densities are in excellent agreement. Below these depths the profiles are parallel but separated by a 2 or 3 m interval. The very high densification rate at Little America V, which begins at 35 m and continues through the pore close-off depth, has been attributed to local areas of high stress Reference Gow(Gow, 1968). Similar sudden departures between observed and predicted profiles may be seen, although not so clearly, at about the 25 m depth in the Roosevelt Island Dome core and at a depth of 20 m in the C-7-3 core: The depths of these intervals of constant density are in general agreement with the depths of anomalies in the velocity-gradient profiles of compressional wave, “C” depths, reported to occur over widespread areas of West Antarctica Reference Robertson and Bentley(Robertson and Bentley, 1975). It was suggested that seismic anomalies, identifiable only in areas of relatively high accumulation rates, represent a transition to an intermediate stage of densification where a different mechanism of metamorphism predominates Reference Robertson and Bentley(Robertson and Bentley, 1975). However, the slopes of the depth –ln [ρ/(ρ _i−ρ)] lines are identical above and below the depths of the anomalies (see Appendix) implying a uniform densification mechanism. There are no intervals of constant density or other anomalies in the profiles of Greenland sites such as Site 2 and South Dome as would be expected if the identification of an additional densification stage were facilitated by greater accumulation rates.

According to the model, an interval of negligible densification might be produced by a long time period of much lower temperature or by a short time period of high accumulation. A high-accumulation event would seem the more likely to occur in Nature. This concept is supported by the fact that within the uncertainties of the dating techniques the depths of the density anomalies may be regarded as synchronous. Stratigraphically determined dates are a.d. 1880 at Old Byrd Station and a.d. 1890 at Little America V Reference Gow(Gow, 1968). Model age calculated from Equations (9) and 11 are a.d. 1883 and a.d. 1885 respectively. Model depths for the a.d. 1885 horizon at other West Antarctic and Ross Ice Shelf locations are 26 m at Roosevelt Island Dome, 20 m at C-7-3, and 17 m at J-9, in good agreement with depths of the density anomalies. Model ages of the seismic “C” depth anomaly, calculated from the data of Reference Robertson and BentleyRobertson and Bentley (1975) for the Little America V-Byrd and Sentinel traverses give a mean date of a.d. 1882 ±34 years. The high degree of scatter is probably due to the combined uncertainties in the accumulation rate and in the “C” depths.

The depth intervals of constant density and deviations between observed and predicted profiles indicate that accumulation rates were greater by a factor of about ten rather than a constant amount at each station. From the accumulation rates and densities at the anomalous zone, ten times the mean annual accumulation would represent the following depth intervals: 3 m at Little America V, 2.3 m at Old Byrd Station, 2.2 m at Roosevelt Island Dome, 1.7 m at C-7-3, and 1.4 m at J-9. These values are in good agreement with the observed intervals, and the concept explains the failure to observe density anomalies or seismic “C” depths at stations with low accumulation rates. Although the intervals of constant density may in fact be merely due to measurement errors, it will be interesting to see if the effects of the Krakatau (1883) eruption are found at the same depths. Known episodes of explosive volcanic activity, including the Krakatau eruption, have been detected by chemical analysis of the Milcent ice core Reference Herron and Langway(Herron and Langway, 1978).

Summary

Empirical rate equations describing the first and second stages of densification have been derived from Greenland and Antarctic ice-core densities. Using these rate equations it is possible to predict depth–density curves and depth-age relationships in polar ice sheets and ice shelves. The only parameters required are the mean annual temperature, annual accumulation rate, and the initial snow density. These can often be determined in the course of a single 10 m deep core drilling operation. Given the mean annual temperature and second stage depth-density data, it is possible to predict the annual accumulation rate to within about 20% of values obtained by other techniques.

Deviations between observed and predicted depth–density profiles may be related to changing climatic conditions. Abnormally high densities at Site 2 are coincident with prominent melt features indicative of unusual summer warmth. Depth intervals of constant density at Old Byrd, Little America V, and possibly at C-7-3, J-9−Roosevelt Island Dome may reflect a synchronous “event” occurring in the 1880’s when large portions of the Antarctic continent received as much as ten times the normal snowfall within a year or two.