Introduction
Paleoclimatic reconstruction is of societal value, and ice cores are prominent in such studies (e.g. National Research Council (US), 2002). Important paleoclimatic variables on ice sheets include the accumulation rate of ice and the mean annual temperature. Each can be estimated in several ways (e.g. Reference PatersonPaterson, 1994; Reference BradleyBradley, 1999), but limitations on existing methods mean it is useful to consider developing new techniques.
Here, we develop an indicator independently suggested by Reference Lipenkov, Duval, Hondoh, Salamatin and BarkovLipenkov and others (1998) and Reference Alley and FitzpatrickAlley and Fitzpatrick (1999), which was originally inspired by the work of Reference GowGow (1968a). As discussed in those sources, Reference Lipenkov, Ryskin and BarkovLipenkov and others (1999, Reference Lipenkov, Salamatin, Duval, Durand, Ohno and Hondoh2005), Reference Lipenkov and HondohLipenkov (2000) and below, the number-density of bubbles in a sample of bubbly ice records the integrated temperature and accumulation rate over the time for that ice to have formed from snow (decades to millennia at different sites). Measured bubble number-density and a firn-densification model can be used to estimate either the firnification temperature or the accumulation rate if the other is known. For simplicity, we develop the application assuming that paleotemperature is known, but the complementary approach is straightforward.
Physical overview
The rate of the sub-freezing transformation of snow to ice is controlled primarily by the temperature and by the overburden pressure, hence the snow accumulation rate (e.g. Reference GowGow, 1968b), with higher temperatures and faster accumulation rates giving faster transformations. The transformation of firn to glacier ice is complete at pore close-off, where the pore spaces between grains are isolated from the free atmosphere to form bubbles (e.g. Reference PatersonPaterson, 1994).
Grains grow during the transformation at a rate that depends primarily on the temperature. (Additional influences from ice flow or firnification deviatoric stresses and impurities are considered in the discussion below, but are generally believed to be minor.) Reference GowGow (1968a) argued from pioneering data that the geometry of the firn/ice at pore close-off is scale-invariant. The bulk density there is nearly constant (Reference Martinerie, Lipenkov, Raynaud, Chappellaz, Barkov and LoriusMartinerie and others, 1994), so bigger grains produce fewer, bigger bubbles. This postulate is supported by analyses of available data performed by Reference Lipenkov, Ryskin and BarkovLikpenkov and others (1999) and in the present study.
Following pore close-off, many processes act to change grain-size, including normal grain growth, grain splitting by polygonization and nucleation and growth of new grains (e.g. Reference AlleyAlley, 1992; Reference Alley, Gow and MeeseAlley and others, 1995; Reference Li and JackaLi and Jacka, 1999). Hence, grain-size in ice quickly loses ‘memory’ of conditions during firnification. In contrast, the number-density of bubbles can retain such a memory for a long time, because gaseous diffusion between bubbles is slow (Reference Ikeda-Fukazawa, Hondoh, Fukumura, Fukazawa and MaeIkeda-Fukazawa and others, 2001), as is collision of bubbles and (usually) splitting of bubbles (Reference WeertmanWeertman, 1968; Reference Alley and FitzpatrickAlley and Fitzpatrick, 1999). The pore space at the time of pore close-off does not immediately consist entirely of spherical bubbles: bubble number-density does increase below the pore close-off depth as cylindrical bubbles are pinched, which can occur over some 50–80 m below pore close-off at cold sites such as Vostok, Dome F and Dome C, Antarctica (personal communication from V. Ya. Lipenkov, 2004); however, if care is taken to measure bubbly ice where bubbles have become predominately spherical, bubble number-density should be largely conserved until the formation of clathrates. The in situ size of bubbles could be used instead of number-density because of the inverse relation between these quantities; however, number-density is more reliable because it is not affected by the relaxation processes during and following core recovery that change bubble size and occasionally produce size-obscuring fractures (e.g. Reference Shoji, Langway, Langway, Oeschger and DansgaardShoji and Langway, 1985, p. 47). (Formation of clathrates raises additional concerns, which we do not address here.)
Firn densification and grain growth are understood relatively well, and can be simulated accurately using empirical models (e.g. Reference GowGow, 1969; Reference Herron and LangwayHerron and Langway, 1980; Reference Alley, Perepezko and BentleyAlley and others, 1986). The only major additional step in developing bubble number-density as a paleoclimatic indicator is better validating and calibrating Gow’s (1968a) conjecture of self-similarity at pore close-off, at least to within the attainable accuracy of the bubble number-density method proposed here. To do so, we identified 16 sites for which bubble number-density, temperature and accumulation rate data were available, and for which conditions of temperature and accumulation rate were likely to have been relatively constant over the time interval when the measured bubbles were forming. We used 15 of the sites for calibration, reserving the 16th to test the complete model.
We used a forward model of firn densification (Reference SpencerSpencer, 2000; Reference Spencer, Alley and CreytsSpencer and others, 2001) with grain growth following Reference GowGow (1969). We inverted for G, the single value of the ratio of bubble number-density to grain number-density at pore close-off, that allowed the most accurate prediction of the observed bubble number-density from observed temperature and accumulation rate at each of the 15 calibration sites. (G could be called the ‘Gow number’.) To assess accuracy, we used the forward model with this single G and with site temperatures to invert for the accumulation rate that best matched the measured bubble number-density at each of the 15 calibration sites and at the 16th validation site. The residual errors from this procedure are encouragingly small, and do not show obvious dependence on site temperature or accumulation rate.
Data
The 15 sites (3 in Greenland and 12 in Antarctica) that were used to calibrate the present model are listed in Table 1 (data from these sites, except those from the Greenland Icecore Project (GRIP) and NorthGRIP sites and Dome C, were also used for model calibration by Reference Lipenkov, Ryskin and BarkovLipenkov and others (1999)). Collectively, the 15 sites have mean annual temperatures and accumulation rates that span a broad range (216–255 K and 22–500 kg m−2 a−1, respectively); additionally, there are published bubble number-densities for each. A 16th site, Dome Fuji, Antarctica, with a mean annual temperature of ~216 K (–57.3°C; Reference Kameda, Azuma, Furukawa, Ageta and TakahashiKameda and others, 1997) and an accumulation rate of 30 kg m−2 a−1 (Reference Dome-FDome-F Ice Core Research Group, 1998) to 32 kg m−2 a−1 (Reference WatanabeWatanabe and others, 1997), provides a test of the model at a temperature slightly colder than any used in the calibration dataset. Bubble number-densities for GRIP, NorthGRIP and Dome Fuji were reduced by 15% to avoid including microbubbles in the calibration (personal communication from V. Ya. Lipenkov, 2004) (microbubbles and their fraction are discussed in Reference Lipenkov and HondohLipenkov (2000)).
Data used for calibration. Different sources sometimes list different temperature, accumulation rate or bubble number-density for a site; where this occurred, we tested combinations of published values as shown. Percent error is the relative difference of the model-derived accumulation rate from the published value
Model
We discuss grain growth and its coupling to a firn-densification model, and then the bubble-number/grain-number ratio, G, at the time of pore close-off.
Grain growth
Reference GowGow (1968a) reported that average grain area in polar firn increases linearly with age
(much as with grain growth in metals (Reference Cole, Feltham and GillamCole and others, 1954)) where 〈r(t)〉 is the average grain-size (m) at time t(a), 〈r 0(t 0)〉 is the average grain-size at time t 0 and k(T) is the crystal growth rate (m2 a−1) at temperature T (K). Reference GowGow (1968a) further recognized that k(T) can be approximated assuming a standard Arrhenius dependence on temperature
where k 0 is a constant (m2 a 1), Eg is the grain-growth activation energy (kJ mol−1) and R is the gas constant (8.31447 kJ mol−1 K−1).
The classic determination of activation energy for grain growth is that of Reference GowGow (1968a). His dataset is plotted in Figure 1. A regression line yields k 0 = 67.4 ± 17.4 m2 a−1 (P <0.05) and Eg = 46.9 ± 4.8 kJ mol−1 (P <0.05), and we use these values. The Reference Chen and SpaepenChen and Spaepen (1991) modification of Equation (2) does not significantly affect the results (see Reference SpencerSpencer, 2005) and so is not adopted here.
Grain-growth activation energy. Data from Reference GowGow (1968a).
We follow the data and conclusions of Reference GowGow (1969) and assume that 〈r 0(t 0)〉2 is a function of temperature, T (K),
We use this extrapolated grain area at time zero, instead of a reasonable approximation to the actual average grain area at the surface, as a simplification. Near the surface, grains grow rapidly in size; however, Reference GowGow (1969) observed that a common grain-size of 0.45′ 10″6 m2 is reached at 3–5 m depth in the firn column for five sites covering a broad range of mean annual temperature and accumulation rate (224–256K and 70–400 kg m−2 a−1, respectively). We also performed a model calibration in which we integrated densification from the surface but grain growth only below a depth of 4 m, at which depth we assumed an average grain area of 0.34 × 10−6 m2 for a set of sites with climates spanning the full range of our dataset, which produced the same results as those obtained using Equation (3) and integrating grain growth from the surface. Certain special sites with extremely low accumulation rates and strong katabatic winds (Reference Fahnestock, Scambos, Shuman, Arthern, Winebrenner and KwokFahnestock and others, 2000) may have anomalously large grains near the surface that are not consistent with our model, but our assumptions are probably quite accurate for most ice-core sites.
Coupled grain-growth/firn-densification model
The firn-densification model used here is that described in Reference Spencer, Alley and CreytsSpencer (2001), where the firn-densification rate equations of Reference Herron and LangwayHerron and Langway (1980) and Pimienta (Reference Barnola, Pimienta, Raynaud and KorotkevichBarnola and others, 1991; Reference Schwander, Sowers, Barnola, Blunier, Fuchs and MalaizéSchwander and others, 1997) were used.
We ran forward models of firn densification with a grain-growth subroutine until the firn pore volume reached its close-off volume, Vc, with the temperature dependence
(Reference Martinerie, Lipenkov, Raynaud, Chappellaz, Barkov and LoriusMartinerie and others, 1994), where Tc is the temperature at the pore close-off depth. This is combined with the weak temperature dependence for the density of ice (Reference BaderBader, 1964; see Reference SpencerSpencer, 2005) to obtain the close-off density.
Bubble-number/grain-number ratio, G
Grain-size at the pore close-off depth, calculated as described above, was converted to grain number-density, N g, assuming spherical grains. (Note that any other assumed shape would yield a slightly different numerical value of G but would not affect the accuracy of the overall calculations.) We then estimated an optimum value of the bubble/grain number-density ratio, G, by minimizing the error between model-implemented accumulation rates and independently estimated accumulation rates, using published bubble-number densities and mean annual temperatures to drive the forward model. We do not present a physical explanation for the value or meaning of G here; we simply make use of the empirical observation of its modern nature and value and postulate that it can be extended to reconstruct paleoclimates.
Results and Discussion
Encouragingly, a single value of G = 2.02 ± 0.08 (P < 0.05) is appropriate for the 15 sites in the calibration dataset, and this value of G works well at the 16th site used for testing. Figure 2 shows the error between published and modeled accumulation rates as a function of mean annual temperature for the 15 sites in the calibration dataset. Because multiple and slightly different values have been published for temperature, accumulation rate and/or bubble number-density for some of the 15 calibration sites, we conducted calculations for the range of published values, giving more than 15 points on Figure 2 (see Table 1). We use the variance of the data plotted in Figure 2 to estimate the uncertainty in the bubble number-density indicator, as the combined uncertainty resulting from estimates of bubble number-density, accumulation rate, temperature, grain growth, grain-size, density and densification rate is otherwise unknown. No trend of G with temperature is evident.
Fractional error of accumulation-rate estimates based on bubble number-density relative to independently derived published values. No functional dependence on temperature is evident. A constant bubble-number/grain-number ratio of 2.02 ± 0.08 (P < 0.05) minimizes the error.
Using G = 2.02 in the forward model with measured site temperature to estimate accumulation rate is accurate to within 41% (P < 0.05) of accumulation-rate estimates derived from independent methods for the set of 15 sites (Fig. 2). Were we forced to pick a single temperature and accumulation rate for each of those sites with multiple published values, we believe that consideration of the timeaveraging lengths and other factors would lead to a set producing a similar value of G but with a smaller error, as described by Reference SpencerSpencer (2005); however, some of the selection criteria would, of necessity, be partially subjective, so here we report the full results even though they are less favorable to the model.
Figure 3 shows both modeled and published annual accumulation rate vs mean annual temperature for the sites used in this study. Also plotted in Figure 3 for reference is the best log–linear fit to the full set of published values of accumulation as a function of temperature. Simply estimating accumulation rate from site temperature and this simple regression is less accurate than estimating using bubble number-density in our model (±71% (P < 0.05) for simple regression vs ±41% (P < 0.05) for our model).
Bubble model evaluation. The modern published accumulation rate is plotted with open circles against published temperature for each site with the accumulation rate for that temperature estimated from our model and measured bubble number-density (filled circles). Some workers have assumed that accumulation depends exponentially on temperature; the best-fit line to the published data is shown for reference. We have not attempted any correction for offsets between cloud and surface temperatures. See Table 1 for data references.
We applied the bubble number-density model to Dome Fuji, Antarctica, a site not part of the calibration dataset and slightly outside its temperature range, with a mean annual temperature of –57.3°C. In Table 2 we show that using the average value for published bubble number-density in Holocene ice from Dome Fuji, plus and minus one standard deviation (neglecting the uppermost measured value because of its proximity to the pore close-off depth), the average result of the present model predicts accumulation rates to within ~6% of the independently derived estimates appearing in the literature.
Test data from Dome Fuji. Percent error is the relative difference of the model-derived accumulation rate from the published value
An additional test is provided by Holocene trends in central Greenland. Reference Pauer, Kipfstuhl, Kuhs and ShojiPauer and others (1999) reported an increase over time in bubble number-densities in the GRIP ice core, from approximately 220 bubbles cm−3 about 5000 years ago (–5 kyr) to 330 bubbles cm−3 recently. Figure 4 is a map of bubble number-densities formed under steady-state climate conditions, plotted as a function of mean annual temperature and annual accumulation rate, with a range of allowed values for GRIP over the most recent 5 kyr indicated. Using the steady-state results of Figure 4 as a guide, the trend in bubble number-density over the most recent 5 kyr in central Greenland is consistent with either an increase in accumulation rate or a decrease in temperature or some combination of the two. One acceptable history would have accumulation and temperature at –5 kyr about 25% less and about 2 K more than recently, respectively. These changes have the same sign as those reconstructed independently (Reference Johnsen, Dahl-Jensen, Dansgaard and GundestrupJohnsen and others, 1995; Reference Cuffey and ClowCuffey and Clow, 1997). The independently reconstructed climatic changes across this interval for the summit of Greenland indicate cooling of 1.5–2 K (Reference Cuffey, Clow, Alley, Stuiver, Waddington and SaltusCuffey and others, 1995; Reference Johnsen, Dahl-Jensen, Dansgaard and GundestrupJohnsen and others, 1995; Reference Dahl-JensenDahl-Jensen and others, 1998; Reference Alley and FitzpatrickAlley and others, 1999) and accumulation-rate increase of ~5% (Reference Cuffey and ClowCuffey and Clow, 1997). The agreement, although not perfect, is encouraging, and falls well within the combined uncertainties. If bubble number-density in the GRIP core records past climate as argued here, an accumulation-rate decrease there over the last 5 kyr would require there to have been a temperature decrease of more than 5 K, well beyond some estimates for the uncertainty in reconstructed temperature change over this time (Reference Dahl-JensenDahl-Jensen and others, 1998), which lends support to the conclusion that there was both an increase in accumulation rate and a decrease in temperature at GRIP over the past 5 kyr.
Bubble number-density reached under steady-state conditions. The triangle between the isolines for 330 and 220 bubbles cm−3 represents possible climatic change over the past 5 kyr allowed by our model and measured bubble number-density. The lower right corner of the triangle represents modern conditions.
Sources of error
Many sources of error are possible in this study. We believe that some issues, including non-conservation of bubble numbers from splitting or coalescence, and failure of the grain-growth or densification models from impurity effects or from excessive vapor transport associated with megadune fields, will be important under certain recognizable situations but not generally. Were this not the case, consistent results on grain growth and bubble number-density, such as in Reference GowGow (1968a, Reference Gow1969), would not have been possible.
The biggest issues are related to the small size of the datasets available. Because the physically based rate equations for firn densification and grain growth are calibrated empirically, quality and quantity of data are critical. Additional datasets beyond the 16 considered here would reduce uncertainties in the model constants. Collection of the available datasets spanned decades and involved different observational techniques, possibly introducing small errors. The datasets also were prepared using different techniques to correct for the offset between feature sizes on section planes and true three-dimensional sizes (Reference Stephenson and ListerStephenson and Lister, 1959; Reference GowGow, 1969). Additional uncertainties may arise from errors in reconstructing climate during the firnification time of the calibration dataset, and from inadequacies in our physical understanding (e.g. failure to include impurity effects on grain growth). Taken together, these raise sufficient questions that our results should not be considered conclusive. Additional work, especially including systematic measurements of physical properties of additional ice cores, would be of great value.
Accumulation rates and grain-sizes are significantly affected by some climatic and topographic conditions, such as megadune fields (Reference Fahnestock, Scambos, Shuman, Arthern, Winebrenner and KwokFahnestock and others, 2000) and areas with glazed surfaces (Reference Frezzotti, Gandolfi, La Marca and UrbiniFrezzotti and others, 2002) caused by hiatuses in accumulation. Clearly these areas fall outside the range of the conditions used to constrain the firn-densification and grain-growth models used here, and would therefore fall outside the applicable range of the current bubble model.
Bubble number-densities used here are assumed to have been conserved above the clathratization zone in all cases. Bubble number-density increases may be brought about by bubble splitting, but such a phenomenon requires large bubbles or high stress, both of which are likely to be encountered only under anomalous conditions (Reference Alley and FitzpatrickAlley and Fitzpatrick, 1999). Bubble number-density reduction through bubble coalescence is considered negligible, as bubble coalescence requires higher strain conditions (Reference WeertmanWeertman, 1968) than are observed above the clathratization zone in all but the most extreme cases.
Grain growth is known to be affected by high quantities of impurities (e.g. Reference Duval and LoriusDuval and Lorius, 1980), and may be affected slightly by lower concentrations (Reference Alley and WoodsAlley and Woods, 1996). However, the effect is typically small (Reference Alley and WoodsAlley and Woods, 1996). Additional consideration may be required in comparing ice age and Holocene samples with different levels of impurities, but we suspect that this effect is almost always small. The sign is known, and the magnitude can be estimated at least crudely (Reference Alley and WoodsAlley and Woods, 1996).
Additional considerations
Any study such as this raises a number of questions. These include issues of clathrates, firn geometry and microbubbles. Conversion of bubbles to clathrates, and existence as clathrate possibly for very long times in deep Antarctic ice, may compromise the bubble memory of firnification processes. Our very preliminary examination gives some hope that the technique can be extended into clathrate-bearing ice (cf. Reference Lipenkov and HondohLipenkov, 2000), based on data including the similarity in trends between bubble number-densities over the most recent deglaciation and clathrate number-densities over the previous deglaciation at Dome Fuji (Reference NaritaNarita and others, 1999). However, issues including clathrate-crystal growth (Uchida and others, 1994; Reference Pauer, Kipfstuhl, Kuhs and ShojiPauer and others, 1999) and clathrate fragmentation (Reference Kipfstuhl, Pauer, Kuhs and ShojiKipfstuhl and others, 2001; personal communication from J. Kipfstuhl, 2003) will require careful consideration before quantitatively accurate and reliable estimates are possible.
Reference Lipenkov and HondohLipenkov (2000) identified a separate distribution of microbubbles in the Vostok core that form in the shallower sections through sublimation–condensation; however, he noted that, above the clathratization zone, microbubbles can be distinguished from the type of bubbles of interest here (those formed by the reduction of pore volume through densification). The extent to which microbubbles affect bubble number-density as a paleoclimatic indicator is unknown at present and will require additional investigation.
We find 2.02 ± 0.08 (P < 0.05) bubbles per grain at pore close-off, but we know of no compelling physical explanation for this value. Bubbles form at four-grain (or greater) boundaries. Several models of space-filling polyhedra have ratios of four-grain boundaries to grains that exceed 1, so we are not surprised to have found a ratio here that is greater than 1. The validity of assuming scale-invariant geometry for firn is lent additional credibility by the similarity in profiles of grain-contact area and coordination number with density in firn columns for different firn types from three sites (Ridge BC and Upstream B, Antarctica, and Site A, Greenland) (Reference AlleyAlley, 1987).
We are investigating the combined interpretation of bubble number-density and of firn thickness, as recorded in gravitational fractionation of trapped gases (Reference Sowers, Bender, Raynaud and KorotkevichSowers and others, 1992). Both are physically based indicators of ice-sheet temperature and accumulation rate averaged over the firnification time. However, the indicators exhibit different dependencies on temperature and accumulation rate, and so produce intersecting (though not orthogonal) isolines of allowable paleoclimatic conditions in temperature/accumulation-rate space. Joint interpretation of firn thickness and bubble number-density thus allows estimation of both paleotemperature and paleo-accumulation rate (albeit with low accuracy), or refinement or validation of other estimates. Additionally, independent paleotemperature estimates combined with paleo-accumulation-rate estimates from both modeled bubble number-density and gas-isotope gravitational fractionation may constrain past convective- and diffusive-zone thickness.
Conclusions
In bubbly glacier ice, where bubble number-density is dynamically stable, there are approximately two bubbles for every grain that existed at the time of pore close-off. Our model, driven by measured bubble number-densities, and estimates of mean annual temperature for modern sites accurately predict independently estimated accumulation rates to within 41% (P < 0.05). Extension of the modern relation between bubble number-density and climate to the last 5 kyr of the GRIP ice-core record is qualitatively consistent with the temperature/accumulation-rate trend estimated with independent methods.
Based on the limited dataset considered here, ice-core bubble number-densities can provide accurate estimates of accumulation rates from temperature histories. Alternatively, bubble number-densities can provide estimates of paleo-temperatures from accumulation-rate histories.
Acknowledgements
We thank all those researchers and support staff who gathered samples, measured grains, counted bubbles or otherwise provided data and insight without which this paper would not have been possible. We also thank J. Kipfstuhl for providing helpful data, and V.Ya. Lipenkov and an anonymous reviewer for many helpful suggestions and corrections. This research was supported in part by the US National Science Foundation Office of Polar Programs through grants including 0087160, 0229609 and 9615554, and by the Comer Foundation.
