Introduction

Knowledge of climate history is of value in understanding and attributing climate change (e.g. Reference SolomonSolomon and others, 2007). Ice cores contribute in important ways, supplementing and extending instrumental data (e.g. Reference Schneider and SteigSchneider and Steig, 2008; Reference Steig, Schneider, Rutherford, Mann, Comiso and ShindellSteig and others, 2009). Data from ice cores are especially valuable because they provide multiple, independent constraints on the history of key climatic variables such as temperature (e.g. Reference Cuffey and PatersonCuffey and Paterson, 2010).

Ice-core techniques for estimating past temperatures include the interpretation of stable-isotope ratios of ice (δ¹⁸O or δD; e.g. Reference Epstein, Sharp and GowEpstein and others, 1970), occurrence of melt layers (e.g. Reference Das and AlleyDas and Alley, 2008), thermal fractionation of gases under temperature gradients across firn, and depth or age of firn (e.g. Reference Severinghaus, Sowers, Brook, Alley and BenderSeveringhaus and others, 1998), together with the inversion of temperatures measured in boreholes (e.g. Reference Cuffey, Alley, Grootes, Bolzan and AnandakrishnanCuffey and others, 1994). All of these have a strong basis in physics but are influenced by processes other than simply temperature change.

Reference Spencer, Alley and FitzpatrickSpencer and others (2006) showed that bubble number-density in polar ice can be quantitatively modeled as a function of temperature and accumulation rate. Here we apply their model to new bubble number-density data obtained from an ice core at the West Antarctic ice sheet (WAIS) Divide site to produce a new paleoclimate reconstruction for the ∼two millennia prior to ∼1700 CE. The WAIS Divide drilling site (Fig. 1) is located at 79°28.058′ S, 112°05.189′ W, ∼160 km from the previous Byrd ice-core drilling site and ∼24 km from the current West Antarctic ice-flow divide (on the Ross Sea side). This drilling location was chosen because it provides an excellent hightime-resolution analogue to the central Greenland ice cores, in terms of ice accumulation rate (∼22 cm a⁻¹), thickness (∼3465 m), average annual surface temperature (−31.1°C), and gas-age–ice-age difference (∼225 years) (Conway and others, unpublished).

Fig. 1. Map of West Antarctica showing the WAIS Divide, Byrd and Siple Dome ice-core drilling locations. The WAIS Divide WDC06A ice-core site is located at 79°28.058′ S, 112°05.189′ W, ∼160 km from the Byrd ice-core drilling site, ∼24 km from the current West Antarctic ice-flow divide and ∼1640 km from McMurdo station. The ice at the WAIS Divide site is ∼3465 m thick, with a bed elevation of ∼1700 m below sea level. Contours represent surface elevation (m). Figure modified from Conway and others (unpublished).

Bubble Number-Density Paleoclimatology

The technique of Reference Spencer, Alley and FitzpatrickSpencer and others (2006) allows reconstruction of either paleotemperatures or paleo-accumulation rates by fitting a semi-empirical steady-state model to measured number-densities of bubbles in (ice-core) glacier ice (see also Reference Lipenkov, Duval, Hondoh, Salamatin and BarkovLipenkov and others, 1998; Reference Alley and FitzpatrickAlley and Fitzpatrick, 1999). Temperature and accumulation rate are the primary drivers of polar firn densification, and the integrated effects of these two drivers on density and grain growth are recorded in the number-density of bubbles formed during the sintering process as the firn is transformed into ice at the pore close-off depth (Reference SpencerSpencer, 1999; Reference Spencer, Alley and FitzpatrickSpencer and others, 2006). That number-density is then conserved in the bubbly ice following pore close-off. This allows use of a combined firn-densification/grain-growth model to invert for either the firnification paleotemperature or paleo-accumulation rate, provided the other parameter is known, as described by Reference Spencer, Alley and FitzpatrickSpencer and others (2006). Here we define pore close-off depth as the depth at which firn has reached a total density of ∼90% that of the average ice density (Reference Herron and LangwayHerron and Langway, 1980).

Accumulation rate can be determined independently for the WAIS Divide site, as described below, so we use bubble data to solve for temperature history here. The transformation of snow to firn, and ultimately to ice, is principally governed by the temperature and by the weight of overlying accumulated snow (Reference GowGow, 1968b). In higher-accumulation environments, this process is achieved more rapidly due to greater overburden pressures. The crystal size at the depth of bubble trapping is controlled by the time to transformation, and by the crystal growth rate, which is primarily controlled by temperature. For these reasons, faster transformation is favored by higher temperature and/or higher accumulation rate. Faster transformation in turn gives a smaller gas-age–ice-age difference, and climate records from bubble number-density with higher time resolution.

Because of this dependence on accumulation rate and temperature, there are large geographic variations in pore close-off depths and firnification times. At higher-accumulation sites such as WAIS Divide, Antarctica, or Summit, Greenland, pore close-off takes only a few hundred years, and may even be achieved in decades at sites of exceptionally high accumulation such as parts of Law Dome, Antarctica. In contrast, at a low-accumulation site such as Vostok, Antarctica, the process can take a few thousand years and is significantly more sensitive to slow temporal accumulation-rate changes (i.e. rate changes through interglacial/glacial transitions) (Reference Sowers, Bender, Raynaud and KorotkevichSowers and others, 1992). The captured bubbles in glacier ice thus preserve a record of environmental conditions during the time that the enclosing ice was still firn, with time resolution that depends on accumulation rate and temperature. Note, however, that the pore close-off process is not instantaneous, introducing slight additional smoothing. For WAIS Divide, pore close-off depth is identified at ∼823 kg m⁻³ ice density (90% of average ice density for the study sample set).

Reference GowGow (1968a) argued that the geometry of firn at pore close-off is essentially self-similar, meaning that larger grains at pore close-off produce fewer but larger bubbles. During burial following pore close-off, the ice grain size itself quickly becomes an unreliable indicator of firnification conditions due to various processes (growth, nucleation and splitting of grains) that depend on the ice deformation rate and not just on temperature and accumulation rate (Reference Alley, Gow and MeeseAlley and others, 1995; Reference Spencer, Alley and FitzpatrickSpencer and others, 2006). However, bubble number-densities in ice reliably preserve information about firnification conditions for a longer time, due to the slow diffusion between bubbles (Reference Ikeda-Fukazawa, Hondoh, Fukumura, Fukazawa and MaeIkeda-Fukazawa and others, 2001, Reference Ikeda-Fukazawa2005) and the near absence of bubble coalescence or splitting in most ice-sheet situations (Reference Alley and FitzpatrickAlley and Fitzpatrick, 1999). It is important to note, however, that bubbles are known to preserve this information only until they collapse into clathrate-hydrate crystals under sufficiently high burial pressures (e.g. Reference Shoji and LangwayShoji and Langway, 1982; Reference Salamatin, Hondoh, Uchida and LipenkovSalamatin and others, 1998; Reference Pauer, Kipfstuhl, Kuhs and ShojiPauer and others, 1999); it remains unclear whether original bubble number-density can be estimated quantitatively from clathrate number-density or from bubbles formed by relaxation of clathrate-containing ice.

For the WAIS Divide site, the bubble–clathrate zone was observed to begin at ∼1100–1200 m depth, which equates to an age of only ∼5 ka because of the high accumulation rate. However, this high accumulation rate allows the technique to give relatively high time resolution, averaging over the ∼225 years to reach the close-off depth at the modern accumulation rate. Here we concentrate on the more recent part of the bubble record, because the older section is in ‘brittle ice’, with air pressures in bubbles sufficiently high to cause fracturing during sample preparation; we plan to extend measurements after the brittle ice has relaxed enough to be analyzed properly.

Methods and Data

In order to employ the bubble number-density model and create the temperature reconstructions in this study, both bubble number-density data and accumulation rates were necessary. Bubble number-density data were acquired by preparing, digitally imaging and analyzing bubble thin sections obtained from ice at 23 depths from the WDC06A WAIS Divide ice core. Accumulation rates necessary for the model were estimated from ice-layer thickness after correcting for ice-flow strain and densification. For dating purposes, a depth–age scale for the WDC06A based on annual layers resolved in chemistry data was used.

Ice-core bubble thin sections for measurement of bubble number-densities were prepared either in the field or at the US National Ice Core Laboratory (NICL)’s −26°C sample-preparation room. First, single-piece 10 cm long samples were cut from the appropriate WDC06A ice-core sections; typically, both vertical and horizontal sections were prepared. The vertical sections cannot be oriented unambiguously relative to the ice-flow direction, but the near isotropy of bubbles observed in horizontal sections means that no bias is introduced by this difficulty when we use a typical bubble size to correct for cut-bubble effects in section analysis, as described below. Samples were collected at ∼20 m depth intervals. In total, the 23 samples from the WDC06A ice core (cut on-site during the 2007/08 WAIS Divide field season) were taken at depths ranging from 120 to 562 m. Before cutting and preparing the bubble thin sections, these 10 cm samples were first used for total ice-density measurements.

Densities were measured by submerging the 10 cm ice samples on a metal ‘basket’ into a small reservoir of isooctane, measuring the volumetric displacement and making the necessary buoyancy corrections to allow for the effect of the basket (Reference GowGow and others, 1997). Measured densities increased smoothly to 915 kg m⁻³ at ∼260 m. The technique is potentially accurate to ±0.3 kg m⁻³, but variability between successive deeper samples was ∼±1 kg m⁻³, perhaps owing to microcracking marking the very earliest influence of brittle-ice behavior. Changes in porosity from bubble compression are estimated to be sufficiently small below 260 m (Reference GowGow, 1968a) that we adopted 915 kg m⁻³ for deeper samples, together with a spline fit for samples between 120 and 260 m. For samples from 0 to 120 m (400–900 kg m⁻³), densities were determined on-site by the volume–mass technique (personal communication from T. Sowers, 2008).

Following density measurement, a vertically oriented slab of ice ∼1.5 cm thick was cut from the side of the sample using a bandsaw, followed by a ∼1.5 cm thick horizontally oriented piece from the top of the sample. Each of these smaller samples was then cleaned using a microtome, sandwiched between glass plates, and bonded using a bead of water on the perimeter of the sample (use of glue on potential bubble sections was avoided because of the possibility of introducing bubbles in the glue).

After setting, each sample was cut in half using a bandsaw to create two separate thin (∼5 mm) samples bound on one side by glass, with one sample destined for thin-section analysis and the other for bubble analysis. The bubble sections were microtomed to a thickness of ∼1.25 mm to ensure that they were thin enough that interpretation would not be complicated by numerous bubbles overlapping one another. Samples thinner than ∼1.25 mm are prone to breakage during microtoming. Finally, sample ice thicknesses were determined by taking measurements of overall sample thickness (including glass plates) and glass-plate thickness alone using a precision digital vernier caliper. Each sample was photographed, with scale, using a Nikon D80 camera (and 18–135 mm Zoom-Nikkor lens) at the NICL on a side-lit stage in a dark room.

Raw images were processed to obtain bubble number-density data (Fig. 2). An area of ∼150 mm² was selected on each sample that was most free of smudges, cracks or other marks and thus provided the cleanest and most consistent area of the bubble section. This size was chosen in order to provide enough bubbles to prevent statistical fluctuations in the results, but without generating too much work for the analysts. Assisted by an edge-detection algorithm, rare remaining non-bubbles (i.e. smudges, glass marks, cracks, etc.) were shaded out and removed, while partial bubbles (those for which the illumination had not produced a completely closed convex figure in the image) were manually edited to elucidate more clearly ambiguous edges, or removed if located on one of the image edges (see Equation (1) below). Next, each edited image was converted into an eight-bit binary image, using the same threshold in all cases, in order to highlight the bubbles and improve identification of bad pixels and overlooked spots. Finally, colored versions of the binary images were overlain on their originally edited counterparts to check for any outstanding errors. Extensive checking shows that our thresholds are set such that we do not introduce spurious ‘bubbles’, so that there is no need to conduct a second check by overlaying the original image on the binary image.

Fig. 2. Four-step ‘clean-up’ process used on ∼150 mm² bubble sections: (a) raw data; (b) non-, partial and edge bubbles blacked out, with ambiguous bubble edges drawn in for clarity (examples noted); (c) after binary conversion, with bubbles highlighted for clarity and bad pixels edited out; and (d) the final step of overlaying the binary on (b) to check for errors.

This complete process was applied to two non-overlapping areas of each raw bubble-section image, and the average was used in our calculations. In addition, the entire image-processing procedure was repeated independently by a different analyst for each of these areas on each image to establish reproducibility. After all final image reads and edits were completed, software scripts were used to automatically measure all of the bubble features and pertinent data. In addition to bubble number-density, various other bubble statistics were measured for potential future use (e.g. average inscribed bubble radius, bubble elongation, bubble perimeter, etc.).

The bubble number-density is the count of bubbles with centers that fall within a specified volume. However, because bubbles are not geometrical points, some bubbles with centers outside a sample volume will be cut by the sides or the upper or lower faces of the sample, and thus will appear within the sample. Corrections must thus be made (e.g. Reference UnderwoodUnderwood, 1970). The corrections are simplified by the fairly narrow size distribution exhibited by the bubbles, allowing us to use the mean bubble size, taken as the average inscribed radius, R, returned by the software. (This introduces a very minor error because the bubbles centered outside the sample influence the average value, but, based on our experience, use of a full Saltykov-type correction would introduce much greater error (Reference UnderwoodUnderwood, 1970; Reference AlleyAlley, 1987).) We eliminated all bubbles touching the sides of our L × L square sampling regions (where L is the edge length), which is equivalent to counting only those bubbles with centers at least a distance R inside the squares. We cannot similarly tell which bubbles have centers above the upper face or below the lower face of a sample, so we count all bubbles with centers in a volume extending a distance R above the upper face and R below the lower face. Using measured ice thicknesses and including these corrections, the final calculated bubble number-density can be represented by

(1)

One additional minor correction was needed. Bubbles form at approximately atmospheric pressure, and are then compressed fairly rapidly until the internal bubble pressure nearly equals the ice overburden pressure. Further compression is then very slow to maintain the internal pressure at the ice overburden pressure, until clathrate formation removes the bubbles (Reference GowGow, 1968a). The finite thickness of a sample, together with the progressive vertical compaction of the ice column from this bubble compression, causes the number of bubbles with centers within a given thickness to increase slightly with depth. This is not a climatic effect and must be corrected for. Here a correction is made to the ice density of 915 kg m⁻³ using

(2)

where B _O is the observed bubble number-density (bubbles cm⁻³), ρ _M is the measured ice density for the given sample (kg m⁻³), ρ _A is the calculated average ice density (kg m⁻³) and B _C is the corrected bubble number-density (bubbles cm⁻³). (In all cases, this correction is quite small. Note also that the incompressible strain of ice from ice-sheet flow increases the horizontal distance between bubble centers while decreasing the vertical distance. Thus the net effect on bubble number-density is zero, so no correction is needed.)

For each horizontal and vertical sample, four separate readings (two readings each for two independent analysts) were averaged together to arrive at a single corrected bubble number-density value; the mean and standard deviation are plotted against depth in Figure 3. The 1σ error in bubble number-density (bubbles cm⁻³) is ±16.78 horizontally and ±13.74 vertically.

Fig. 3. Final horizontal and vertical bubble number-densities plotted against depth. Error bars are 1σ uncertainty of reproducibility across samples.

Computations and Results

As discussed by Reference Spencer, Alley and FitzpatrickSpencer and others (2006), the bubble number-density depends on the integrated effect of the accumulation rate and temperature on firn densification and grain growth over the time interval between snow deposition and bubble close-off. An ultimate goal is to develop fully time-dependent inverse techniques for paleoclimatic reconstructions. However, when climate changes are relatively small and slow (as indicated post facto for the WAIS Divide site, and as expected for most sites during the Holocene), a useful short cut is to apply forward modeling to map the dependence of bubble number-density on constant accumulation rate and temperature spanning all possible combinations with degree and millimeter precision. Then an average accumulation rate over a specific interval can be used in conjunction with this map to estimate the average temperature over the time from deposition to bubble formation. Pending development of a time-dependent inverse technique, we follow the pseudo-steady-state approach here, using the accumulation–temperature–bubble mapping of Reference Spencer, Alley and FitzpatrickSpencer and others (2006). We further estimate accumulation rates averaged over 57 m (ice) thickness, the approximate modern value for the depth to pore close-off, assumed constant over time, and corrected for thinning from ice flow. In times of strong accumulation-rate change, this assumption could cause use of a significantly erroneous average value, but post facto it introduces little error.

Accumulation rate can be estimated from ice-flow correction of an accurate depth–age scale; at the shallow depths considered here, the uncertainties in the flow correction are small (e.g. Reference Cuffey and PatersonCuffey and Paterson, 2010). Annual layers in the WAIS Divide ice core have been counted using various techniques including visible identification, dielectric properties analysis and chemical analysis (e.g. Reference TaylorTaylor and others, 2004; Reference Banta, McConnell, Frey, Bales and TaylorBanta and others, 2008; Reference McGwire, McConnell, Alley, Banta, Hargreaves and TaylorMcGwire and others, 2008). The techniques generally agree well. The chemistry-based depth–age scale is believed to provide the most accurate representation of annual-layer timing and thicknesses, and was therefore used in this study.

Raw annual-layer thicknesses from the chemistry-based depth–age scale were corrected first for density, and second for plastic deformation by way of an ice-flow model with constant strain rate with depth (Reference NyeNye, 1963):

(3)

where T_ρ is the density-adjusted layer thickness (cm), H is the total ice-sheet thickness at WAIS Divide (3465 m), Z is the height of the given layer from the bed of the ice sheet (m) and T _O is the final ice-flow and density-corrected layer thickness (cm), or ice equivalent accumulation rate. Layer thicknesses (and thus ice accumulation rates) are plotted against their age over the last ∼2400 years (Fig. 4). At greater depths, a more accurate flow model will be required, but little error is introduced here. The general trend toward lower accumulation rate with time agrees with the results of Reference Neumann, Conway, Waddington, Catania and MorseNeumann and others (2008). The greatest potential source of error is in the dating, but we believe that errors are <1%. Uncertainties in accumulation due to density, thickness and ice-flow correction were too small to appear clearly in figures or notably affect the calculations.

Fig. 4. WAIS Divide accumulation rates calculated after density and ice-flow strain correction. The black curve is a 225 year running average.

For each horizontal and vertical sample in this study, a paleotemperature was determined using the bubble number-density, average accumulation rate, and modeled steady-state values of Reference Spencer, Alley and FitzpatrickSpencer and others (2006; see their fig. 4 for graphical representation of the full tabular output). Little change in bubble number-density occurs across the dataset (Fig. 3). The trend toward lower accumulation rate with time (Fig. 4) favors fewer, larger bubbles. The constant bubble number-density thus indicates a synchronous cooling, which favors more but smaller bubbles. This result is evident in the calculated temperatures as shown in Figure 5. The model error for determination of absolute temperature from Reference Spencer, Alley and FitzpatrickSpencer and others (2006) is ±0.9°C. Our counting reproducibility, and thus our measurement error contributing to uncertainty in temperature changes, is ±0.5°C. Regression of the results shown in Figure 5 yields trends over the complete dataset of ∼1.65°C and 1.67°C for horizontal and vertical samples, respectively, with 95% uncertainty of ∼1°C. There is no reason to believe the trend is linear, but inspection of the data, or comparison of points from early and late in the history, indicate a decrease of slightly less than 2°C. For now, we take our regression result as our best estimate, indicating a decrease of ∼1.66 ± 1°C from ∼257 BCE to ∼1686 CE.

Fig. 5. Reconstructed temperatures plotted against age. Both horizontal and vertical samples show a similar decreasing linear trend. Error bars for temperature are derived from the bubble number-density standard deviation errors. ‘Error bars’ for age represent the averaging time for each sample. The linear regressions yield an average cooling of ∼1.668°C.

For comparison with the bubble number-density temperature estimates, we analyzed the ¹⁸O/¹⁶O ratios of ice in 50 cm long contiguous samples of the WDC06A core. The measurements were done at the University of Washington using equilibration at 35°C with CO₂ (e.g. Reference CraigCraig, 1961) with a Micromass Isoprime mass spectrometer with ‘Aquaprep’ sample preparation system. Values are referenced to working standards calibrated to the International Atomic Energy Agency standards GISP (Greenland Ice Sheet Precipitation), SLAP (Standard Light Antarctic Precipitation) and VSMOW (Vienna Standard Mean Ocean Water), and reported using standard δ¹⁸O notation, expressed as a deviation (‰) from VSMOW. External accuracy of the measurements is <0.1‰.

The δ¹⁸O data are shown in Figure 6, and are converted to approximate temperatures using the current average annual surface temperature for WAIS Divide of ∼−31.1°C, an approximation for the modern δ¹⁸O average value of −33.6‰ and an assumed δ¹⁸O proxy calibration value of ±0.8‰ representing : ±1°C (Reference MorganMorgan, 1982). Agreement between the isotopic and bubble-number-density temperature trends is excellent.

Fig. 6. Running-average temperature histories from δ¹⁸O (225 year average) and bubble number-density (six-sample). In each case, the temperature error bar is the 1σ uncertainty of the average, and the horizontal bar is the averaging length. Shaded bands are added for ease of viewing, and to identify possibly interesting trends.

Discussion

Nearly constant measured bubble number-density is observed over time (Fig. 3). At constant temperature, the steady decrease in accumulation rate over time (Fig. 4) would have caused longer firnification times, hence formation of larger but fewer bubbles (Reference SpencerSpencer, 1999). A decrease in temperature, slowing grain growth and thus reducing bubble size, increasing number-density and counteracting the effect of decreased accumulation, is thus indicated (Reference GowGow, 1969).

A decrease in accumulation rate often results from a decrease in temperature through dependence on the saturation vapor pressure, at ∼7%°C⁻¹ (Reference Denton, Alley, Comer and BroeckerDenton and others, 2005; Reference Banta, McConnell, Frey, Bales and TaylorBanta and others, 2008). The relation observed in the results from WAIS Divide is 9%°C⁻¹ (Fig. 7), consistent within the combined uncertainties.

Fig. 7. Reconstructed average temperatures plotted against measured accumulation rates. Plot shows a ∼9% increase in accumulation rate per °C warming, broadly consistent with control by saturation vapor pressure. Error bars for temperature are derived from the bubble number-density standard deviation errors.

The reconstructed temperature and accumulation rate applied where and when the ice was deposited. However, because both climate and the ice sheet can change, and ice is carried by flow even if the ice sheet is not changing, a reconstructed change in temperature or accumulation rate may have multiple causes. In our case, consideration of the modern flow field indicates that the reconstructed changes primarily reflect climate change. We first consider the nonlocal origin of the deeper ice analyzed, and then possible effects of changes in surface elevation over time.

The WAIS Divide site is ∼24 km on the Ross Sea side of the ice divide, so modern flow is bringing ice to the core site that was deposited as snowfall uphill and closer to the Amundsen Sea. Reference Conway and RasmussenConway and Rasmussen (2009; see also Reference Neumann, Conway, Waddington, Catania and MorseNeumann and others, 2008) reported the modern flow velocity at the site as ∼3 m a⁻¹, decreasing up-glacier. Thus 2400 year old ice has come from 6–7 km up-glacier where the modern surface elevation is ∼5 m higher than at the drill site. Surface slopes in the vicinity of the core site show some variability (Reference Conway and RasmussenConway and Rasmussen, 2009), and today the slope near the core site is not as steep as in regions both closer to the ice divide and down-flow of the core site. However, calculations using these steeper nearby slopes and the modern flow velocity yield a deposition site for 2400 year old ice only ∼10 m higher than at the core site. Assuming a lapse rate of 1°C (100 m)⁻¹ or slightly less, modern ice flow will have introduced an apparent warming of <0.1°C over the length of our climate record, small compared to the observed signal. The modern accumulation rate decreases along ice flow in response to regional gradients, so the accumulation for the oldest ice in our sample was almost 5% larger than for the youngest ice (Reference Conway and RasmussenConway and Rasmussen, 2009). Thus, about one-third of the apparent accumulation-rate change over time may have arisen from flow through the spatial accumulation-rate pattern, with the other two-thirds representing climate change.

At present, the site of the WAIS Divide core is nearly in balance, neither thickening nor thinning (Reference Conway and RasmussenConway and Rasmussen, 2009). However, 20 km down-flow toward the Ross Sea the ice is thickening by 0.1 m a⁻¹, and 30 km up-flow toward the Amundsen Sea the ice is thinning by 0.1 m a⁻¹ .If a rate of 0.1 m a⁻¹ through a lapse rate of 1°C (100 m)⁻¹ had applied to the WAIS Divide site over the entire history of our climate record, then, depending on the sign, it could eliminate the entire climatic signal or double it. But, with the modern behavior at WAIS Divide, almost none of our signal was caused by changing surface elevation over time.

The ice divide is currently migrating toward the Ross Sea at ∼10 m a⁻¹, faster than the 3 m a⁻¹ flow of the ice (Reference Conway and RasmussenConway and Rasmussen, 2009). If the entire thickness-change pattern has migrated with the divide, then the WAIS Divide site experienced elevation gain slightly less than 0.1 m a⁻¹ at the beginning of our record, dropping to zero recently. The resulting elevation change through the lapse rate might account for half of our reconstructed cooling. Given the large uncertainties, we have not attempted to use a detailed flow model to correct for ice flow, but we do not find a straightforward explanation of the entire signal based on ice flow.

An alternate approach is to examine the behavior of the site in a simulation of the deglaciation and Holocene in a whole-ice-sheet model. Recent three-dimensional West Antarctic ice sheet modeling results (Reference Pollard and DeContoPollard and DeConto, 2009; personal communication from D. Pollard, 2009) indicate surface-elevation changes of <100 m over the time interval corresponding to our study. Different model runs initialized in slightly different ways produce distinct trends, with both thickening and thinning observed, but consistently remaining small. Reference SteigSteig and others (2001) also found only small changes at the ice divide for a variety of model initial conditions. Strong conclusions should not be drawn from these modeling results, but they support our extrapolation from the data of Reference Conway and RasmussenConway and Rasmussen (2009) that at least half and probably most of the reconstructed cooling represents a climatic change rather than an ice-flow effect.

Our reconstructed cooling has the appropriate sign for the expected effects of decreasing radiative equilibrium temperature in response to the effect of decreasing obliquity on summer duration (Reference Huybers and DentonHuybers and Denton, 2008). There is no necessary conflict with the melt-layer results from relatively nearby Siple Dome showing midsummer warming over the same interval (Reference Das and AlleyDas and Alley, 2005, Reference Das and Alley2008), because those likely reflect midsummer conditions that are especially sensitive to midsummer insolation intensity, whereas the bubble number-density reflects mean annual temperature that is more sensitive to obliquity control on the seasonal distribution of intensity (Reference Huybers and DentonHuybers and Denton, 2008).

Notably, however, δ¹⁸O data from both the nearby Byrd (Reference Johnsen, Dansgaard, Clausen and LangwayJohnsen and others, 1972) and Siple Dome (Reference SchillerSchiller, 2007) ice-coring sites appear to indicate a slight warming trend over the same period as in this study. Reference SteigSteig and others (2001) argued that Holocene trends at Byrd reflect ice-sheet thickness changes moving the surface through the atmospheric lapse rate, as well as the effects of horizontal ice flow causing older samples to have been deposited at higher, colder sites. Changes of order 100 m would be sufficient to significantly shift surface-temperature trends from climatic trends. Independent evidence points toward larger elevation changes at Byrd and Siple Dome than at WAIS Divide (Ackert and others, 1999, Reference Ackert, Mukhopadhyay, Parizek and Borns2007), such that our interpretation of a primarily climatic cooling over the study interval at WAIS Divide does not necessarily indicate a different trend from those other sites.

Conclusions

The new paleoclimatic indicator based on ice-core bubble number-density developed by Reference Spencer, Alley and FitzpatrickSpencer and others (2006), when applied to samples taken from the WAIS Divide WDC06A ice core, and with annual-layer thicknesses calculated from a chemistry-based depth–age scale, shows an approximately linear cooling of ∼1.7°C over the ∼two millennia prior to ∼1700 CE. This cooling is consistent with δ¹⁸O isotope data for the WDC06A ice core. Additional work is required, but, either by extrapolating the modern ice-flow pattern back in time or by modeling the response of the sampling region to the deglaciation and Holocene, we find that a significant part of the reconstructed cooling represents climate change rather than the effects of ice flow. In turn, this cooling is consistent with Milankovitch obliquity forcing. The drop in temperature at WAIS Divide was accompanied by a decrease in accumulation rate, broadly consistent with dependence on saturation vapor pressure.