Introduction

Understanding the densification of polar firn is integral to several areas of glaciological research. For example, mass-balance estimates for ice sheets derived by monitoring surface elevation (e.g.Reference Wingham, Ridout, Scharroo, Arthern and ShumWingham and others, 1998) may be affected by a contribution arising from the response of the firn column to climate change. The thermal properties of firn are highly density-dependent (Reference YenYen, 1981), so the reconstruction of a temperature history based on direct borehole measurements (paleothermometry; e.g. Reference JohnsenJohnsen, 1977) is facilitated by knowledge of the firn column’s depth–density history. Knowledge of the pore close-off depth (where the air becomes trapped in the network of pores and is isolated from the atmosphere above) is required to calculate the age offset between the ice and the gas it contains (e.g. Reference Schwander, Sowers, Barnola, Blunier, Fuchs and MalaizeSchwander and others, 1997). For these reasons and others, several empirical or semi-empirical firn-densification models have been developed (e.g. Reference Herron and LangwayHerron and Langway, 1980; Reference AlleyAlley, 1987a; Reference Barnola, Pimienta, Raynaud and KorotkevichBarnola and others, 1991; Reference Arnaud, Lipenkov, Barnola, Gay and DuvalArnaud and others, 1998). A range of such models is valuable, to help quantify model-dependent uncertainties, and because different models may perform better in different temperature, accumulation, density or other ranges.

Based on previous models and on hot-isostatic-pressing experience, we assume that densification rates are functionally dependent on temperature, density and overburden load, and that power-law creep is the dominant mechanism. Parameters in such a model are difficult to determine accurately from first principles owing to geometric and physical uncertainties. Instead, we leave all the parameters free to vary during a data inversion using 38 depth–density profiles from locations collectively spanning a large range in mean annual temperature and accumulation rate (216–256 K and 0.022–1.2 m w.e. a⁻¹, respectively). The set of best-fit rate equations contains non-physical values for some parameters and is seen to fit poorly near the surface; however, the rate equations appear more consistent with the data at greater depths than some current dynamic models appearing in the literature that use a higher-order constitutive form (e.g. Reference Barnola, Pimienta, Raynaud and KorotkevichBarnola and others 1991; Reference Schwander, Sowers, Barnola, Blunier, Fuchs and MalaizeSchwander and others, 1997), particularly towards the extremes of mean annual temperature and accumulation rate of the dataset.

Methods

We searched the literature to find complete datasets, including all of depth (z), density (ρ(z)), accumulation rate and mean annual temperature (T _10m). This yielded a total of 92 sites (Reference SpencerSpencer, 2000). Inspection of the data raised two issues. First, some of the datasets had potential problems. These include sparse depth sampling, relatively large quoted errors, surprisingly large scatter of densities in comparison to the majority of profiles, and smoothing before publication without details of how this was accomplished. Second, many stations had almost identical temperature and accumulation rate, such as the numerous cores collected during site selection for and operation of the GRIP and GISP2 sites in central Greenland. Use of all the data in such clusters would tend to bias any empirical model to fit those conditions extremely accurately at the expense of fitting the wide range of temperatures and accumulation rates actually observed on ice sheets. To solve these problems, as described in Reference SpencerSpencer (2000), we selected a 38-site subset of the data to obtain high data quality and coverage across a broad range of and T _10m. The locations and climate characteristics of the 38 sites selected to constrain empirical rate equations are presented in Table 1. (We did conduct our analyses on the whole dataset as well as on the 38-site subset, with broadly similar results, but the full dataset yielded a model with poorer ability to simulate depth–density profiles, especially at the extremes of temperature and accumulation rate such as at Vostok, Antarctica.)

Table 1. Site characteristics for the 38 firn depth–density profiles used to constrain Equation (2)

Depth–density profiles were smoothed using spline methods (Reference SpencerSpencer, 2000). Sampling intervals on some of the depth–density profiles included in our 38-site subset were rather irregular, and smoothing served to regularize these. The smoothed profiles were converted to yield densification rates by making the usual assumption that each site’s observed depth–density profile has reached steady state with the observed mean annual temperature (generally taken as the 10 m firn temperature, T _10m) and observed accumulation rate; thus, , where ρ is the density (kg m⁻³), t is time (years), is the accumulation rate (kg m⁻² a⁻¹) and z is the depth from the surface (m).

The error associated with the assumption that the firn profiles used are in a steady state with the observed climate is difficult to assess. It will vary from site to site because the time for a firn column’s development can range from dozens to thousands of years depending on site characteristics. In general, Holocene climates have been relatively stable compared to glacial–interglacial changes (e.g. Reference LoriusLorius and others, 1985), and the time for firn to change to ice is short compared to the duration of the Holocene at all sites, so errors associated with the steady-state assumption are likely to be small. These assumptions do not preclude that the sites used have experienced a secular trend large enough to affect the method followed here; however, without the steady-state assumption, the depth–density profiles can no longer be viewed as representing a trajectory for parcels of firn, and inverting the data to constrain the rate equations directly would be impossible. Instead, one might use an iterative approach whereby a known climate record is applied to a forward model, the density profile of which would be compared with site data at each iteration; this method is left for future investigations.

We adopt the standard method of dividing the densification rates into three stages, with decreasing densification rate and increasing density with increasing depth. The transition from the first, lowest-density stage to the second is associated with the achievement of random closest packing of grains (Reference Anderson, Benson and KingeryAnderson and Benson, 1963) and thus the near-cessation of grain rearrangement by boundary sliding (Reference AlleyAlley, 1987a). This transition density is typically taken to be 550 kg m⁻³ (Reference BensonBenson, 1962; Reference Herron and LangwayHerron and Langway, 1980) and we follow this here. The transition from the second to the third stage occurs where the pore space becomes isolated from the atmosphere above and the increase with densification in the air pressure in the isolated pores opposes further densification. By combining empirical relations for the pore volume found in firn at the isolation depth (Reference Martinerie, Lipenkov, Raynaud, Chappellaz, Barkov and LoriusMartinerie and others, 1994) and the density of pure ice (Reference Schwander, Sowers, Barnola, Blunier, Fuchs and MalaizeSchwander and others, 1997), a temperature-dependent close-off density of (kg m⁻³) was used to define the transition from the second to the third stage, where T _c is the temperature (K) at the pore close-off density and typically today is quite close to T _10m.

Densification rates were assumed to have a standard Arrhenius dependence on temperature (dρ/dt ∝ exp(−E/RT)). Here E is an activation energy (kJ mole⁻¹), R is the gas constant (kJ mole⁻¹ K⁻¹) and T is the temperature (K).

The densification rate in the upper few meters is increased by the large temperature variability there. This is because the exponential dependence of densification rate on temperature means that the summertime increase in densification rate above that appropriate for the mean annual temperature is larger than the wintertime decrease (i.e. the average of an Arrhenius term experienced over a full cycle of a time-varying temperature is greater in magnitude than the Arrhenius term calculated using the average temperature). This effect is only significant in the upper few meters, because of the rapid attenuation of temperature variability with depth in firn.

We corrected for this effect by increasing the site temperatures in the upper few meters to an effective tempera ture, T _eff, obtained by integrating the Arrhenius term over a sinusoidally varying annual cycle:

(1)

(Reference SpencerSpencer, 2000), assuming an activation energy of 50 kJ mole⁻¹ based on available data from other studies. A(z) is the attenuated temperature amplitude with depth, A(z) = A ₀ exp(−0.457z), where A ₀ is the surface temperature amplitude. We obtained this A(z) in two ways, which agreed to two significant figures. In the first, we forward-modeled the near-surface firn using the Reference Herron and LangwayHerron and Langway (1980) first-stage rate equation and the empirical relation from men (1981) for the temperature- and density-dependent thermal conductivity of snow using mean conditions of temperature and accumulation rate for our dataset. In the second, we took observed depth-density profiles spanning our dataset and modeled the temperature with depth for a specified surface amplitude (Reference PatersonPaterson, 1994). Both approaches will produce small errors for sites with very high accumulation rates (e.g. DE08 on Law Dome, with an accumulation rate of 1.2 m w.e. a⁻¹ (see Table 1) and a vertical firn velocity of 2.5 m a⁻¹, which is high enough to slightly increase the downward velocity of the annual temperature cycle and thus increase the densification rate to a few meters depth); we make no additional correction for this, which is not significant for any of our sites other than DE08 on Law Dome. Figure 1 shows how much the effective temperature is increased above the mean annual temperature as a function of the activation energy assumed, the temperature amplitude and the mean annual temperature. The effect is significant in the upper ∼3 m but not below, as shown.

Fig. 1. Effective temperature increase over the mean annual temperature vs activation energy for mean annual temperatures of 200–250 K and seasonal, sinusoidally varying temperature amplitudes of 5–20 K. Activation energies of 50 kJ mole⁻¹ were used in inversions. Taking 20°C as the surface temperature amplitude, 15°C amplitude corresponds to about 0.6 m depth, 10°C to 1.5 m and 5°C to 3 m. Clearly, the effect is not large deeper.

The data with effective temperatures for 38 sites were then used to constrain the rate equations

(2)

where ρ is the firn density (kg m⁻³), ρ _ice is the density of ice (here kept constant at 917 kg m⁻³ for simplicity), P _eff is the effective pressure (Pa), and C _1–5 are the 15 constants (5 constants × 3 stages) determined by inverting the data using a generalized version of Newton’s method for those that can be made linear and a simplex method for those that cannot. The effective pressure, P _eff, was assumed to be proportional to the overburden pressure and inversely proportional to both the relative density (ρ/ρ _ice) and the relative grain-contact area (fraction of grain surface in contact with ice rather than air). Correction for relative density is required because a given load supported on pillars of ice causes a higher stress than if supported on a continuous slab of ice; the additional correction is required because narrow necks between grains in these ice columns further amplify the load. The relative grain-contact area as a function of density was assumed to be equal to the best-fit, second-order polynomial of the stereographically determined grain-contact area measurements presented in Reference AlleyAlley (1987b). Throughout the third stage the relative density and the average grain-contact area both approach 1, eliminating any effective-pressure increase over that of the overburden alone; also, bubble pressure (over that of the air pressure at the time of pore close-off) was subtracted from the overburden pressure for densities above the pore close-off density.

The functional form of Equation (2) was derived by Reference Wilkinson and AshbyWilkinson and Ashby (1975) for hot-isostatic power-law creep around both spherical pores and cylindrical channels. Hot-isostatic power-law creep-rate equations have been applied by others to firn (Reference Barnola, Pimienta, Raynaud and KorotkevichBarnola and others, 1991; Reference Arnaud, Lipenkov, Barnola, Gay and DuvalArnaud and others, 1998), but with decreasing usefulness at lower density. There are two main reasons for the poorer fit at lower densities: the contribution of other densification mechanisms (e.g. grain-boundary diffusion or grain-boundary sliding) and the increasing deviation from a state of hydrostatic stress. However, Equation (2) was chosen for its simplicity and lack of any variables not readily accessible through paleoclimatic reconstruction (such as grain or neck size, wind speed, etc.). Also, grain-boundary sliding, which likely dominates densification near the surface (Reference AlleyAlley, 1987a), should increase with temperature, load and inverse density as in Equation (2). In Wilkinson and Ashby’s derivation of Equation (2), the constants C _3–5 are equivalent, but here they were left free to vary, to be constrained only by the data.

Results and Discussion

The best-fit coefficients for Equation (2) are displayed in Table 2. While C ₁ and C ₂ have the appropriate sign and magnitude, the value of C ₅ is negative throughout stages 1 and 2, being positive only during stage 3 but still much lower than the value of 1–3 that might be expected (Reference PatersonPaterson, 1994) from the Reference Wilkinson and AshbyWilkinson and Ashby (1975) model. Negative C ₅ indicates that an increase in pressure causes a decrease in densification rate. This clearly non-physical result is a consequence of decreasing rates of densification with increasing depth (and therefore decreasing rates of densification with increasing pressure), and it demonstrates that the density portion of the constitutive form is counteracting the non-physical exponent of the pressure term. Pressure and density both increase with depth, so they are not independent variables. Any inversion to simultaneously constrain the density and pressure components of the constitutive form will automatically include the best-fit density dependence for the pressure term, as well as the best fit to the true constitutive form that is available with the constitutive form actually being used. Additionally, inadequacy of the assumption that each depth–density plot represents a trajectory for a parcel of firn (i.e. that the system is in steady state) could account for a significant deviation of the coefficients from their expected values. For example, consistent deviation of past accumulation rates from those of today (perhaps because of a Little Ice Age signal) amongst enough sites would result in the coefficients being different from idealized expectations. However, we believe that the lack of independence between depth and pressure is more important in controlling the results of the inversion.

Table 2. The best-fit coefficients of Equation (2) as determined by the 38 sites described in Table 1

Equation (2) with the best-fit coefficients appearing in Table 2 is plotted with the data in Figure 2 for eight sites that collectively span the full range of mean annual temperature and accumulation rate of the dataset. We consider the Pimienta model (Reference Barnola, Pimienta, Raynaud and KorotkevichBarnola and others, 1991) to be the current state of the art in firn-densification modeling (e.g. Reference Schwander, Sowers, Barnola, Blunier, Fuchs and MalaizeSchwander and others, 1997), and for comparison it is also plotted in Figure 2. The present model consistently predicts monotonically decreasing rates of densification with increasing density, in accord with observations and in contrast to some previous models in certain circumstances. The present model overestimates the densification rate compared to both the data and the Pimienta model at Byrd, Antarctica (using the data of Reference GowGow (1968) from “old” Byrd Station), and underestimates the densification rate at Law Dome, Antarctica. However, the present model is seen to match the data well at sites slightly warmer and cooler than both Byrd and Law Dome.

Fig. 2. Densification-rate density profiles for eight sites: (a) Vostok, Antarctica; (b) Dome C, Antarctica; (c) South Pole, Antarctica; (d) Dome GRIP, Greenland; (e) Byrd, Antarctica; (f) Camp Century, Greenland; (g) Law Dome, Antarctica; and (h) Dye 3, Greenland. The densification-rate data were calculated from the depth–density data and accumulation rate at each site (see Table 1 for references to data sources). The Pimienta model (Reference Barnola, Pimienta, Raynaud and KorotkevichBarnola and others, 1991) is plotted with the present model (Equation (2) with coefficients from Table 2) for comparison.

The poorer fit for the present model at Law Dome is to some extent a consequence of the anomalously high accumulation rate experienced there. As can be seen in Figure 3, the accumulation rate for Law Dome is well above the Arrhenius temperature-dependence trend. This in turn means that Law Dome has anomalously high loads at a given density compared to most polar sites. The inversion matches the majority of sites, and so is not especially accurate for Law Dome. The Reference Herron and LangwayHerron and Langway (1980) model that explicitly includes accumulation rate matches Law Dome better. However, load is physically more directly relevant to densification than is accumulation rate and so might be expected to perform better overall in formulating a model. Indeed, although our model does not match Law Dome as well as it might, overall our model’s performance is somewhat better than that of Reference Herron and LangwayHerron and Langway (1980) with its accumulation-rate dependence.

Fig. 3. Accumulation-rate Arrhenius temperature plot for 86 sites. The filled circles are the 38 sites used in the data inversion. With an R ² value of 0.48 for the best-fit, first-order line, the Arrhenius dependence accounts for slightly less than half of the variation in accumulation rates. Law Dome is seen to plot significantly above the best-fit line.

For modern Vostok climate (mean annual temperature of 216 K and accumulation rate of 22 mm w.e. a⁻¹) the present model predicts that a firn density of 830 kg m⁻³ is reached at 94 m depth, equal to the observed depth, whereas the Pimienta model (using the Herron and Langway first-stage rate equation) predicts the same density is reached at 105 m depth. However, when extrapolating to mean annual temperatures below the range found in the dataset, the present model predicts a significantly deeper firn column than does the Pimienta model. Using a mean annual temperature of 204 K and an accumulation rate of 10 mm w.e. a⁻¹ to simulate the Vostok climate at the Last Glacial Maximum, the Pimienta model predicts that a density of 830 kg m⁻³ is reached at 137 m depth, whereas the present model predicts that the same density is reached at 181 m. Simple considerations using the Reference Wilkinson and AshbyWilkinson and Ashby (1975) model and assuming that strain rate increases with the cube of the stress at high stress magnitudes suggest that firn is unlikely to delay close-off of bubbles to such great depth.

The greater firn-column heights predicted by the present model using mean annual temperatures outside the range of the dataset substantially eliminate confidence in extrapolations using the present model. We believe the disparity between the two models at lower temperatures is due partly to the use of the Reference Herron and LangwayHerron and Langway (1980) firststage rate equation in the Pimienta model of Reference Barnola, Pimienta, Raynaud and KorotkevichBarnola and others (1991) (which has an activation energy smaller than the present model’s first-stage activation energy by a factor of approximately 4). However, as revealed in Figure 2, the Pimienta model already begins to reveal non-physical trends at both the high and low extremes of mean annual temperature in the dataset, producing sometimes erratic jumps in densification rate that do not match observations.

Conclusions

Because we have used independent techniques and partially independent data, the differences between the Pimienta model and ours provide some information on uncertainty in firn-densification modeling. We believe that our model typically does a slightly better job for conditions spanned by our input dataset, at the cost of producing much poorer results for extrapolation beyond our input conditions. The greater firn-column heights predicted by the present model using mean annual temperatures lower than the range of the dataset, and the non-physical trend of the model used for comparison near the temperature extremes of the dataset, lessen confidence in firn-column heights predicted for climates outside the range of modern values.