Large Arctic Temperature Change at the Wisconsin-Holocene Glacial Transition

Analysis of borehole temperature and Greenland Ice Sheet Project II ice-core isotopic composition reveals that the warming from average glacial conditions to the Holocene in central Greenland was large, approximately 15°C. This is at least three times the coincident temperature change in the tropics and mid-latitudes. The coldest periods of the last glacial were probably 21°C colder than at present over the Greenland ice sheet.

The Greenland Ice Sheet Project II (GISP2) deep ice core has yielded a remarkable history of the oxygen isotopic composition of central Greenland snowfall [8180 of ice (1)] that extends through the last glacial period (2). The nearby Greenland Ice Core Project (GRIP) record (3) is essentially identical for ice formed after the 110,000-year-old Eemian interglacial, and both are similar to isotope histories obtained in other Greenland cores, giving confidence that these cores record aspects of regional climate (4). Using both empirical data (5) and physical models for isotope fractionation (6), paleoclimatologists have interpreted 8180 to be a measure of environmental temperature T at the core site, through a simple relation that we call the isotopic paleothermometer: 8QO = aT + 1, where ao and L are constants. There are two obstacles to making this interpretation sound. First, the coefficients a and 1 are not known a priori (7-9) because many factors in addition to local environmental temperature affect isotopic composition. These include changes in sea-surface composition and temperature (10), changes in atmospheric circulation ( 11), changes in cloud temperature, which may be different from changes in surface temperature (12), changes in the seasonality of precipitation (13), and postdepositional isotopic exchange in the snowpack (14). Second, all of these factors may vary through time in such a way that a single, linear relation between strong motivation to seek paleotemperature information that is entirely independent of isotopic history (15,16) to calibrate the paleothermometer. We have obtained such information by measuring temperature at depth in the ice sheet, and we use this information to evaluate ao and 1.
During the summer of 1994, one of us (G.D.C.) measured temperature in the 3044m-deep GISP2 core hole from 70 m below the surface to the base of the ice sheet. At that time, the thermal perturbation from drilling had decayed to less than 0.04°C, so the temperature in the borehole matched the temperature in the surrounding ice sheet at this accuracy and better (17). To determine the coefficients ax and 1 in the isotopic paleothermometer, we used the GISP2 8180 record and an initial guess for (a and 13 to specify a 100,000-year history of environmental temperature. We then calculated subsurface temperatures using T as the forcing function on the upper surface of the ice sheet in a linked heat-and ice-flow model. Finally, we adjusted a, 1, and the geothermal heat flux from the underlying bedrock, using the Levenberg-Marquardt method (18) to minimize the mismatch between modeled and measured subsurface temperatures. For this purpose we defined the mismatch index J as a weighted integral over ice depth z of the squared difference between modeled and measured subsurface temperatures (M and 0, respectively) Here UD and (rt are weighting functions that assign relative importance to misfit in various parts of the borehole: a, assigns more weight to the upper part of the borehole, where the temperature does not depend on poorly known ice dynamical quantities, and a, assigns more weight to the lower part of the borehole, where the temperature depends on longer intervals of the surface temperature history. The parameter f, which is adjustable, controls the trade-off between these opposing weighting schemes SCIENCE * VOL. 270 * 20 OCTOBER 1995 (19). The solution that minimizes J is nonunique because we are free to choose f.
The heat-flow component of our model is a numerical solution to the advectiondiffusion equation with heat sources (20). The model is one-dimensional (vertical), with a movable upper boundary to allow changes in the ice sheet thickness. These changes, and the vertical ice velocity responsible for heat advection, are calculated with standard glaciologic assumptions for flow on the flank of an ice divide with simple parameterizations to account for two-dimensional effects. In our model, the ice sheet responds to local changes in snow accumulation rate, surface temperature, and ice crystal fabric and to distant changes in ice margin position (21). To calculate vertical heat advection most confidently, and to account for two-dimensional effects on ice particle paths, we tuned the vertical velocity so that the modem depth-age scale matched our model within a small tolerance. Because of this tuning, our conclusions are insensitive to poorly known aspects of the ice dynamics. The GISP2 depth-age scale was determined by annual layer counts to about 40,000 years ago (40 ka) and by correlation to the ocean-core time scale (SPECMAP) through analysis of the oxygen isotopic composition of 02 gas for older ice (22-24).
The snow accumulation rate history b(t) exerts a dominant control on ice sheet thickness and vertical ice flow, and hence vertical heat advection. For the most recent 35,000 years, we derived b(t) from the layer thickness measurements of Meese and others (23), corrected for ice-flow thinning. Our correction uses strain calculations as done previously (24, 25), but we also included the dependence of layer thinning on the thickness history of the ice sheet (26, 27), which in turn depends on the temperature history. Before 35 ka in the model l3 -34 kwt's.  Fig. 1. The central Greenland 8180 history for the most recent 40,000 years. The smooth curve results when this history is filtered to mimic the thermal averaging in the ice sheet (45). All temperature histories that give this same curve when filtered are indistinguishable to borehole thermometry (29). The right axis shows our calibrated temperature scale. runs, we calculate b(t) from the oxygen isotopes using a linear correlation between 8180 and our b(t) for the most recent 35,000 years (28).
Because heat diffusion damps high-frequency temperature changes as they propagate from the surface down into the ice sheet (thermal averaging), information on rapid environmental temperature changes in the past is poorly retained by the presentday temperature field O(z) in the ice sheet (29). Thermal averaging is more extensive for older climatic events. In contrast, the GISP2 518Q record retains information about rapid climatic changes. If we degrade the 8180 record so that it retains only the age-dependent low-frequency content that can be recovered from the present-day temperature field (Fig. 1), then the abrupt termination of the Younger Dryas, the Younger Dryas itself, and the B0lling/Aller0d period become minor features of the history, and earlier interstadial events are no longer evident. Thus, our isotope calibration is sensitive mainly to the long warming from full glacial conditions to the Holocene, and to Holocene temperature changes (30).
We find the optimal linear paleothermometer to be 518Q = 0.327T -24.8 if we assume ice dynamics are well known (f = 1; refer to Eq. 1), and 818Q = 0.335T -24.5 if we assume ice dynamics are poorly known (f = 1000) (31). Using these calibration constants, we find a remarkably good fit between the temperatures measured in the borehole (Fig. 2) and the corresponding modeled temperatures; the model accounts for 99.88% of the variance of the measured profile relative to steady-state. This is strong evidence that 8180 is indeed a faithful proxy for long-term average temperature at this site. There is no better explanation for the success of such a simple calibration, given the small number of free parameters in the inversion.
However, the fit is not perfect. For in-    3) and insensitive to changes in ice dynamical parameters ( Table 1). The average temperature difference between the Wisconsin Glaciation and the Holocene is therefore large (Fig. 1), 140 to 16°C, and the coldest periods of the last glacial were probably 21°C colder than at present (Fig. 1). The climatic deglacial temperature change (at constant elevation) may be 10 to 2°C larger than this because the Greenland ice sheet was probably thinner during the glacial as a result of a substantial reduction in accumulation rate. Geologic evidence suggests that the margins of the ice sheet retracted by about 100 km during the Wisconsin-Holocene transition (35). Using this value, and assuming a symmetrical retreat of east and west margins, we estimate that the ice sheet thickened by 250 m from the last glacial maximum to the present, and at least 100 m from average glacial conditions to the present (36). Recent estimates of the Wisconsin-to-Holocene warming in the low mid-latitudes are 40 to 6°C. This result is based on a variety of methods, including snow line depression studies (37), palynology (38), noble gas paleothermometry applied to ground water (39), and stable-isotope paleothermometry applied to coral reefs (40, 41). The -8°C temperature change commonly inferred from ice-core isotopic records (37), including those from the new GISP2 and GRIP cores (2, 4), using the modem spatial value for a of 0.60 to 0.67 per mil 'C-1, is only slightly larger than recent estimates from the tropics. By contrast, we have shown that the temperature change in central Greenland was three to four times larger than that in the tropics, a result that is consistent with borehole temperature analyses at Dye 3 in southem Greenland (42). Many models have suggested that initially minor changes in global temperature will be magnified in the Arctic, with possibly major consequences for sea level and planetary albedo (43). Our data not only confirm that such amplification happened in the past but also show this amplification to be larger than generally thought. vironmental conditions at GISP2, the system sensitivity is 0.00014°C; the accuracy of the measurements traced to the National Institute of Standards and Technology is -0.0045°C. Temperature measurements were acquired every 2 s while lowering a custom sensor with 20 thermistor beads down the hole at -6 cm s-t. Instrumental noise during this experiment was limited to less than 0.001°C. Data from the moving sensor were then deconvolved to find the stationary borehole temperatures by the methods described in R. W. Saltus and G. D. Clow, U.S. Geol. Surv. Open-File Rep. 94-254 (1994). The precision of the completely processed data is estimated to be better than 0.0010C for wavelengths greater than about 6 m. Drilling activities between 1990 and 1993 disturbed the temperatures in the ice surrounding the borehole. Comparison of temperature measurements acquired in 1994 and 1995 in the deep borehole, and temperatures acquired to 163 m in a nearby air-filled hole, shows that the magnitude of the disturbance in 1994 was about 0.040C near the top of the fluid-filled section of the main borehole. This disturbance diminished to about 0.0010 to 0.0020C at 400 m in depth and remained at this level to the bottom of the hole. Because the disturbance is largely restricted to the upper part of the borehole, its effect on our inferred deglacial temperature change is negligible (about 1.5% of the temperature change rheology, ancient snowfall rates, ice sheet geometry changes, and initial temperature constitute the major sources of uncertainty in calculating subsurface temperatures. This uncertainty is large deep in the borehole and small near the surface, suggesting that we should ascribe more importance to misfit in the upper part of the borehole: crD is a rigorous means of doing so. However, temperatures high in the borehole depend on only the most recent several centuries of the temperature history, whereas temperatures deep in the ice sheet depend on many tens of thousands of years of the temperature history. This suggests we should ascribe more importance to misfit in the lower part of the borehole, to find the isotopic paleothermometer that applies best to the whole isotopic his-  (1995)] also found relatively high values of a for recent times. This trend in a could reflect a changing climate system (32); it could also indicate that 8180 is more sensitive to temperature changes over years to centuries than to temperature changes over periods of millennia or longer. The Cuffey et al. (7) a value is slightly higher than ours (0.53 versus 0.46 per mil 'C-1) because their measurements extend 40 m closer to the ice sheet surface and hence are weighted more heavily to recent temperature changes [see figure 3 in (7)]. 34. Changes in seasonality of snow accumulation may contribute to the difference between the spatial and temporal gradients [(13) and P. J. Fawcett, A. M. Agustsdottir, R. B. Alley, Eos 76 (spring suppl.), 177 (1995)]. We used raw isotopic data from the ice core, with no adjustment for changes in source-water isotopic composition, whereas Jouzel et al. (8) assumed that source-water changes equaled the effect of ice-volume changes on mean oceanic composition, but removed these source-effects from their calculated last glacial maximum 8180 when they estimated the spatial gradient a. Recalculation of the Jouzel et al. (8) temporal gradient with the snowfall isotopic change uncorrected for source effects would yield an even closer match to our result.  (2), (ii) the uncertainties in some key input parameters such as the melting history of ice sheets (3) in the analysis of the postglacial rebound, or (iii) the density-tovelocity conversion factor (4) in the analysis of the geoid. One strategy to get around these difficulties is to combine seismological observations of anisotropy and laboratory studies of deformation-induced lattice preferred orientation [for example (5)]. The anisotropic structure of deformed materials depends on deformation mechanisms (and deformation geometry) [for a review, see (6)], and anisotropic structures can be observed seismologically as far down as the center (that is, inner core) of the Earth (7). Plastic deformation by diffusion or superplastic creep will result in an isotropic structure, whereas deformation by dislocation creep or twinning results in an anisotropic structure. Thus, although this approach will not provide direct estimates of the absolute values of viscosities, it provides information as to the rheological constitutive relation (stress or grain-size dependence of viscosity) and hence indirectly indicates rheological discontinuities or weakening associated with SCIENCE * VOL. 270 * 20 OCTOBER 1995 grain-size reduction (5,8).
Here we apply this strategy to the lower mantle. One of the most striking observations of the lower mantle is the absence of significant seismic anisotropy (9) even though the dominant mineral in the lower mantle, orthorhombic (Mg,Fe)SiO3 perovskite, has significant elastic anisotropy (10,11 ). The absence of observed anisotropy in the lower mantle could be attributed to (i) chaotic convection; (ii) limited plastic deformation; (iii) the anisotropic structure in the lower mantle which happens to be such that, for seismic waves traveling nearly vertically [such as those used in seismological studies (9)], the amount of shear wave splitting is small (12); or (iv) the fact that the deformation does not result in an anisotropic structure. The first hypothesis means that the lower mantle materials could have anisotropic structure (possibly due to deformation by dislocation creep), but the scale of coherent deformation is much smaller than the length of typical seismic wave paths (for ScS or SKS) so that, on average, no appreciable anisotropy would be detected. This hypothesis is unlikely because seismic tomography indicates that the lower mantle structure is dominated by long wave length features (13), and a high viscosity of the lower mantle will make chaotic convection difficult to achieve (14). The second hypothesis is also untenable because seismic tomography indicates the presence of downgoing and upwelling currents in the lower mantle (15), and the Rayleigh number for the lower mantle is likely to exceed the critical value [see, for example, (1)]. Discrimination of the last two alternatives requires an investigation of the relation between the nature of deformation and seismic anisotropy in lower mantle materials.
The most direct data on this subject must ultimately come from high-pressure, high-temperature deformation experiments performed on polycrystalline (Mg,Fe)SiO3