Significance The magnitude and timing of Antarctic temperature change through the last deglaciation reveal key aspects of Earth’s climate system. Prior attempts to reconstruct this history relied on isotopic indicators without absolute calibration. To overcome this limitation, we combined isotopic data with measurements of in situ temperatures along a 3.4-km-deep borehole. Deglacial warming in Antarctica was two to three times larger than the contemporaneous global temperature change, quantifying the extent to which feedback processes amplify global changes in polar regions, a key prediction of climate models. Warming progressed earlier in Antarctica than in the Northern Hemisphere but coincident with glacier recession in southern mountain ranges, a manifestation of changing oceanic heat transport, insolation, and atmospheric CO 2 that can further test models.

Abstract The most recent glacial to interglacial transition constitutes a remarkable natural experiment for learning how Earth’s climate responds to various forcings, including a rise in atmospheric CO 2 . This transition has left a direct thermal remnant in the polar ice sheets, where the exceptional purity and continual accumulation of ice permit analyses not possible in other settings. For Antarctica, the deglacial warming has previously been constrained only by the water isotopic composition in ice cores, without an absolute thermometric assessment of the isotopes’ sensitivity to temperature. To overcome this limitation, we measured temperatures in a deep borehole and analyzed them together with ice-core data to reconstruct the surface temperature history of West Antarctica. The deglacial warming was 11.3 ± 1.8 ∘ C, approximately two to three times the global average, in agreement with theoretical expectations for Antarctic amplification of planetary temperature changes. Consistent with evidence from glacier retreat in Southern Hemisphere mountain ranges, the Antarctic warming was mostly completed by 15 kyBP, several millennia earlier than in the Northern Hemisphere. These results constrain the role of variable oceanic heat transport between hemispheres during deglaciation and quantitatively bound the direct influence of global climate forcings on Antarctic temperature. Although climate models perform well on average in this context, some recent syntheses of deglacial climate history have underestimated Antarctic warming and the models with lowest sensitivity can be discounted.

From the Last Glacial Maximum (LGM) around 21 ka to the middle of the Holocene, increased greenhouse gas concentrations and reduced reflectivities of the surface and atmosphere directly increased the uptake of energy to Earth’s climate system by about 7 W ⋅ m − 2 (1⇓–3) and warmed the surface by 3–6 C ∘ on average (4⇓⇓–7). Although contributing little to this global average because of the comparatively small area involved, the warming in polar regions holds particular interest. In addition to driving changes of ice sheets, permafrost, and hydrology and modulating oceanic and atmospheric circulations, polar warming partly controlled both the evolution of surface reflectivity and the transfer of carbon dioxide from ocean to atmosphere and hence the climate forcing itself. Further, reconstructions of polar warming during deglaciation permit quantification of one key prediction of climate theory—that feedback processes amplify temperature changes in polar regions relative to the global average (4, 8, 9), a phenomenon referred to as polar amplification. Arctic data reveal a warming three to four times the global average based on a wide variety of indicators (6), including combined analyses of ice-core data and borehole temperatures (10, 11). Limited available constraints suggest a smaller but still amplified Antarctic warming, roughly 1.5 to 2.5 times greater than the global average (6, 12). This Antarctic estimate, however, derives from the isotopic composition of ice measured in cores. The scaling between ice isotopic composition and temperature depends on a great many factors (section 15.5 in ref. 13), such as the seasonal timing of snowfall and rate-dependent fractionations in clouds, that remain poorly known and are expected to vary with time. Moreover, the low accumulation rates at East Antarctic core sites have precluded convincing quantifications of the isotopic thermometers using independent information from borehole temperatures, the most direct legacy of past climate. One Antarctic study used such information (14), but the small accumulation rate at the site meant that diffusive heat transport greatly dominated advection, rendering the method imprecise.

Here we present a reconstruction of West Antarctic temperature history from the site of the West Antarctic Ice Sheet Divide ice core (WDC) (15, 16). This site is uniquely suitable for analyses using borehole temperatures, given the combination of thick ice and high snow accumulation rate, together with the wealth of information provided by the 68-ky core record. Elsewhere, similarly favorable conditions prevail only in central Greenland (10).

SI Reconstruction Strategy This section summarizes the reconstruction strategy following the presentation in the main text, but with added details and comments. To determine the surface temperature history T s ( t ) , we optimized the match between temperatures measured in the borehole and those calculated with a numerical model of heat transfer driven by various T s ( t ) scenarios as a boundary condition. To start, we filtered the deuterium ice-isotopic record ( δ D ) to remove high-frequency variability. Using the time derivative of this filtered history ( δ ˙ ), we then optimized the temperature variation (relative to modern) given by Δ T s ( δ ) = ∫ t γ − 1 ( t ) δ ˙ ( t ) d t , [S1]where the coefficient γ ( t ) takes three values (Table S1), one for each of the three major periods of isotopic change (deglacial, early to mid-Holocene, and late Holocene). For the most recent period 0–0.2 ka we do not calibrate γ but instead add and adjust a trend whose shape follows the temperature reconstruction of Orsi et al. (52) derived from analysis of a neighboring shallow borehole. Fig. 2A (main text) shows the result. This calibrated Δ T s ( δ ) is used in all subsequent stages of the analysis, without further adjustments. Next we introduced three basis functions g i ( t ) with ranges ∈ ( 0,1 ) , so that T s ( t ) = T o + Δ T s ( δ ) + ∑ i = 1 3 c i ⋅ g i ( t ) , [S2]and optimized the coefficients c i and the constant T o , the modern temperature. The basis functions reflect the broadest pattern of variations seen in the observed temperature–depth profile and allow for adjustments to average temperatures of the late glacial period, early Holocene, and mid- to late Holocene (Fig. S2). If Δ T s ( δ ) matched perfectly the true temperature history, then we would here find c i = 0 . In fact, nonzero c i are found to improve the model performance significantly (Table S2), indicating that the δ D of ice does not always preserve a linear signal of the long-term evolution of surface temperatures. On the other hand, an optimization of the basis functions without isotopes performs significantly worse than the isotopes alone (Table S2, top entry); the isotopes do record important information about past temperatures. Uncertainty in how the measured temperature matches ice temperature at a given depth is approximately 0.5 cK, but additional unknowns in the model raise the tolerance to 0.8 cK (SI Calculation of Limits and Tolerance section). Any models with rms mismatch differing by less than this amount should be regarded as equivalent, whereas models performing worse by this amount or better compared with the final optimized reconstruction are rejected. Fig. 2A (main text) also shows the history derived from basis functions without isotopic inputs, along with an optimized version prescribing a cooling before the LGM. Note that it would be a mistake to regard these nonisotopic models as better representations of the borehole temperature information. The borehole temperatures do discriminate between different histories containing a wide range of frequency contents, and these nonisotopic versions can be rejected, as can Eq. S1. The “long-term averages” of climatic temperatures preserved in the borehole temperature signal are millennia in the Holocene and tens of millennia in the glacial, but uneven weighting in the averages and the simultaneous constraint provided by the measurements at different depths increase the discriminatory power. The δ 15 N gas-isotope data provide a test of the T s ( t ) estimated using Eq. S2 and an opportunity for improvements. The history of accumulation rate ( b ˙ ( t ) ) can be derived in two independent ways, one from T s ( t ) and δ 15 N , using a model for the evolving firn density and thickness, and a second one from the ice core’s observed annual layer thicknesses, corrected for strain using parameters (total ice thickness, basal melt rate) determined in the optimization of Eq. S2. Discrepancies between the two versions of b ˙ ( t ) can in principle be eliminated by adjusting the temperature history by an amount Δ T N determined by densification physics. Thus, we write T s ( t ) = T o + Δ T s ( δ ) + ∑ i = 1 3 c i ⋅ g i ( t ) + ω ( t ) ⋅ Δ T N ( t ) , [S3]where ω ( t ) ranges between 0 and 1. This coefficient safeguards against possible problems with the method; the accuracy of Δ T N as a thermometer is not well established, especially given imperfections of firn densification models and irregularities of the gas-trapping process at the firn base. We examined reconstructions with various ω values and used borehole temperatures to identify the range of admissible scenarios. Including the term Δ T N in Eq. S3 with ω > 0.5 for all but the late Holocene reduces the rms error of the optimized model temperature profile to 1.47 cK (Table S2). This is an excellent match considering the total range of the observed temperature profile (Fig. 1, main text), but still larger than the 0.8-cK tolerance defined by uncertainties in the measurement and glaciological variables. To assess how much of the remaining mismatch could be eliminated without altering the broad patterns implied by δ D and δ 15 N , we introduced a random perturbation term to Eq. S3 before optimization. The magnitude of perturbations was restricted to one-half the difference between the optimized Eq. S3 scenario and its multimillennial average. Of 60 random scenarios we tested, the best one reduced the rms error to 1.11 cK. This improvement (rms error reduced by 0.36 cK) is not large enough to negate the model without perturbations and regardless is achieved without any significant changes in the reconstructed temperature: less than 0.8 °C warmer during deglaciation and less than 0.05 °C warmer at the LGM (Fig. 2A, main text). Fig. 1 (main text) illustrates the match between model and observed temperature–depth profiles for this best-performing case. Regardless of whether the final perturbations are included, the reconstructed temperature history provides an excellent match to the constraints provided by measured borehole temperatures, ice-core layer thicknesses, and δ 15 N isotopes. Table S2 summarizes metrics of model performance for all stages of the analysis.

SI Ice-Isotopic Data The δ D of ice was measured by laser spectroscopy at the University of Washington as described in Steig et al. (24), and the associated δ 18 O data have been published previously (15, 16, 24). Use of δ 18 O rather than δ D in our calculations makes no discernible difference because δ D scales linearly with δ 18 O . The data are reported as per mille deviations from Vienna Standard Mean Ocean Water (VSMOW), normalized to the VSMOW-Standard Light Antarctic Precipitation standard water scale. The precision of the δ D measurements is better than 0.8%.

SI Nitrogen-Isotopic Data The δ 15 N of N 2 trapped in the ice was measured at Scripps Institution of Oceanography, University of California. Air was extracted from ∼ 12-g ice samples, using a melt–refreeze technique, and collected in stainless steel tubes at liquid-He temperature. δ 15 N was analyzed using conventional dual-inlet isotope ratio mass spectrometry (IRMS) on a Thermo Finnigan Delta V mass spectrometer. Results were normalized to La Jolla (CA) air, and routine analytical corrections were applied (55, 56). The δ 15 N data and model fits were reported previously by Buizert et al. (47) and are shown in Fig. S1.

SI Calculation of Temperatures vs. Depth We use the control volume method (49) to solve the energy balance equation on a grid of 200 nodes in the ice and 25 nodes in subjacent bedrock. Grid spacing is smaller near the surface and increases with depth. The grid is rewritten by interpolation to allow for changes of ice thickness. The Kelvin temperature T evolves with time t and elevation above bed z according to ρ c ∂ T ∂ t + ρ c w ∂ T ∂ z = ∂ ∂ z [ k ∂ T ∂ z ] + S ˙ , [S4]where ρ is density, c is the temperature-dependent heat capacity, w is the vertical velocity, k is the temperature-dependent thermal conductivity, and S ˙ are the source terms. The latter are calculated from the rate of work of ice deformation and firn compaction and are of minor significance. Equations for thermal parameters are those of Yen’s compilation (57), but with thermal conductivity multiplied by a factor a k so that k = a k k Yen . We use a k = 1.02 to conform with recent measurements (58), but vary it from 1.00 to 1.04 in sensitivity tests. Advection is of primary importance because ice accumulates, deforms, and flows downward. The age vs. depth relation for the core (SI Age–Depth Relation and Layer Thicknesses section) provides a tight constraint on the net vertical displacement of layers, making calculated temperatures only weakly sensitive to uncertain details of our flow model. The vertical velocity w (negative in the ice sheet because of downward flow) depends on basal melt rate m ˙ and vertical normal strain rate ϵ ˙ z z , w = − m ˙ + ∫ 0 z ϵ ˙ z z d z [S5] − w ( H , t ) ⋅ ρ ( z ) ρ i = b ˙ − d H d t [S6]for ice-equivalent accumulation rate b ˙ ( t ) and ice thickness H ( t ) . The history of H can take any specified form (see below). The history of b ˙ derives from the ice core’s measured annual layer thicknesses λ ( z ) , b ˙ ( t ) = λ ( z ) ⋅ exp ( − ∫ t ϵ ˙ z z ( z λ , t ) d t ) , [S7]where z λ indicates the depth following the layer through time. With the addition of a correction for densification near the ice surface (the density as a function of depth below the surface is assumed invariant over time), the vertical strain rate follows a Dansgaard–Johnsen model: uniform in the upper part of the ice sheet and varying linearly below a kink height specified by parameter ξ = 0.2 (varied in sensitivity tests), ϵ ˙ z z = ϵ ˙ ∘ + w ρ d ρ d z for ξ H < z < H [S8] ϵ ˙ z z = ϵ ˙ ∘ [ f b + [ 1 − f b ] z ξ H ] for 0 < z < ξ H . [S9] The requirement that Eqs. S5 and S6 match at the surface fixes the value for ϵ ˙ ∘ . Parameter f b defines how the strain rate at the bed compares to the strain rate in the upper 1 − ξ fraction of the ice thickness. f b is sometimes referred to as the “sliding parameter” but this label misleads; the fraction of ice motion due to sliding can change with distance along a flowline, and such a gradient allows f b to take values less than or greater than unity. In our optimized models we use f b = 1.3 , because this value allows b ˙ inferred from ice strain to match that inferred from N 15 in the deepest layers (ages in excess of 45 ka). A value f b > 1 implies that the fraction of motion due to sliding increases along the flowline. Whether this is true or not upstream of the WAIS Divide is unknown. Fortunately, the parameter f b has little influence on our temperature reconstruction, because the optimization adjusts basal melt rate m ˙ . Using a larger value for f b is largely compensated in the temperature profile throughout most of the ice thickness by a smaller value for m ˙ . Reducing f b from 1.3 to 0.7, for example, causes the optimized LGM temperature to decrease by only 0.08 °C. As an initial condition, we set the profile T ( z ) at 120 ka, the starting time of model calculations, equal to the modern observed temperature–depth profile.

SI Basis Functions The three basis functions ( g i ( t ) in Eq. 2, main text) are simple, are piecewise linear, and range from zero to unity (Fig. S2). The g 1 and g 2 center on 4 ka and 10 ka, respectively, the approximate age equivalents of the depths of the two prominent extrema in the borehole temperature profile (Fig. 1, main text). The g 3 extends through the LGM. In the standard case, g 1 is defined by vertex values (0, 1, 1, 0) at ages of (0, 3.5, 4.5, 8) ka, whereas for calculations of limits the vertex ages range as (0–2.5, 2–5.5, 4–6, 6–9) ka. Ranges for vertex ages of g 2 are (4–8, 8–12, 12– ∞ , 15– ∞ ) ka. For g 3 , the ages of the (0, 1) vertices range as (12–15, 15– ∞ ) ka.

SI Firn Thickness Model and Use of Nitrogen Isotope Data Calculation of δ 15 ( t ) . Given simultaneous histories b ˙ ( t ) and T s ( t ) , we calculated the firn density profile ρ ( H − z , t ) from the empirical model of Herron and Langway (20), recast to use overburden load as a driving variable as in refs. 28 and 47. Equations of the densification models are given in appendix A of ref. 47. The density profile, in turn, is used to calculate δ N 15 ( t ) of gas at the depth of trapping by combining the barometric equation (section 15.3 of ref. 13 and ref. 18), an estimate of thermal fractionation due to temperature gradients (19), and a small convective-zone thickness offset. The calculated δ 15 N for current climate matches the observed modern value at the WAIS Divide. Calculation of Δ T N ( t ) . Starting with b ˙ ( t ) derived from observed ice-core layer thicknesses and the strain model derived with the preliminary optimized T s ( t ) , we calculated δ N 15 ( t ) , using the model summarized in the previous paragraph. Comparison with measured δ N 15 ( t ) then defined an adjustment to b ˙ ( t ) , by calculating the derivative of mismatch between model and measured δ N 15 ( t ) with respect to perturbations of b ˙ ( t ) and shooting for a perfect match by linear extrapolation. This process was iterated until measured and modeled δ N 15 ( t ) agreed as well as possible in terms of rms error, thus defining a new adjusted accumulation rate history b ˙ 1 ( t ) . The temperature-history correction Δ T N ( t ) was then calculated as a function of the ratio of the adjusted accumulation rate b ˙ 1 ( t ) to the original accumulation rate b ˙ o ( t ) from strain-corrected layer thicknesses. This function describes the coupled dependence of δ 15 N on temperature and accumulation rate according to firn density models and the barometric equation and has been calibrated against modern-day observations at various sites (figure S1 of ref. 28). We used a simple approximation of the relationship, appropriate for the range of temperatures and accumulation rates at the WAIS Divide, Δ T N ( t ) = ν ln [ b ˙ 1 ( t ) b ˙ o ( t ) ] , [s10]in which the constant ν = 8.475 for temperature change in Kelvins. The entire process was iterated several times to achieve an optimal match against all constraints. Fig. S1A compares the optimized model’s two histories of accumulation rate, one estimated from ice strain and the other from reconstructed temperatures via δ 15 N . We also repeated analyses with two alternative firn densification models (21, 54) and found that they performed not as well as Herron and Langway’s with respect to the borehole temperature match after optimization. The versions that produced an acceptable fit to the borehole temperatures corresponded to surface temperature histories well within the range of uncertainty defined using the Herron–Langway model. Details about use of the alternative densification models are provided in ref. 47. Note that for most of the optimizations in this study, the density profile was assumed to be in equilibrium with the instantaneous values of T s and b ˙ (a quasi-steady state). Only for the final optimizations of Eq. 2 (main text) with and without perturbations, and for the three different densification models, was the densification model run in fully time-dependent mode. The Coefficient 𝝎 . The following method was used to define ω in Eq. 2 of the main text for the nominal optimized models. Considering only the portion of the record with climate similar to modern (12 ka and younger), we perturbed T s in discrete intervals to assess whether, upon optimization, model performance improves or degrades by adjusting in the direction indicated by Δ T N . Improvement was found for all but the latest Holocene, so we took ω = 1 for 3.5–12 ka but ω = 0 for < 3 ka. (In the late Holocene, the reconstruction is already tightly constrained by the combination of ice isotopes and borehole temperatures.) All possibilities in the 0–1 range are included subsequently in tests to define limits on the temperature history. After iterative optimizations of Eq. 2 (main text) and experiments with ω for the > 12 -ka period, we found that in mid-Holocene there remains a difference of ∼ 0.015 m ⋅ y − 1 between the two reconstructions of b ˙ (Fig. S1A), implying that the δ 15 N thermometer does not entirely agree with the borehole temperatures. We assume this is a deficiency of the former and therefore do not force the match between b ˙ versions to be better than 0.015 m ⋅ y − 1 in the pre-Holocene. (There is no reason to expect it to perform better in LGM conditions than in the Holocene when conditions are similar to the present, whereas the cost of allowing a better match of b ˙ is a greater deviation of T s from the ice isotopes. The borehole temperatures average over too great a length at LGM to choose.) This means that ω decreases to ∼ 0.5 for the LGM and prior in our standard case.

SI Calculation of Limits and Tolerance Sensitivity tests reveal seven principal sources of uncertainty with nontrivial effects on reconstructed LGM temperature (ranked in descending order of importance): the ice thickness history, the depth–age scale, the temperature dependence of thermal conductivity of ice, the basal temperature (equivalent to the local melting point of ice), the depth of the borehole temperature sensor, and the residual thermal disturbance from drilling. Setting each of these to maximum and minimum plausible values and reoptimizing define ± 2 σ limits on reconstructed temperatures (Table S4). Adding all of these individual contributions in quadrature then defines the ± 2 σ range due to uncertain inputs. In this exercise it was also found that choosing the right combination of values for these uncertain quantities could improve the rms borehole temperature mismatch by ∼ 0.3 cK. In combination with the ∼ 0.5-cK uncertainty of the direct temperature measurement, this indicates that any reconstructed temperature history that provides an rms mismatch within ∼ 0.8 cK of the mismatch for the standard optimal model should be regarded as equally viable. We produced a set of alternative optimized reconstructions by four processes: introducing random perturbations to Eq. 2 (main text) as described previously; varying the duration and timing of the basis functions g i ( t ) in Eq. S3; varying the coefficient ω ( t ) of Δ T N in Eq. S3; and defining T s as weighted averages η T s 1 + ( 1 − η ) T s 2 , where T s 1 and T s 2 are, respectively, the standard optimal model (Eq. S3) and the calibrated ice-isotope model (Eq. S1), and η ranges from 0 to 1. From this entire set of optimized reconstructions, we accepted as viable those that met three criteria: (i) borehole temperature rms mismatch no worse than 0.8 cK compared with the standard model; (ii) accumulation rates calculated using nitrogen isotopes and strain match over the period 0–30 ka to within an rms difference of 1.36 cm ⋅ y − 1 , twice the rms difference for the standard model; and (iii) model N 15 values, averaged in 1-ka windows, never trend outside the range of measured N 15 points in the Holocene. Considered together, at every age the viable reconstructions span a range of temperatures. We identify the maximum and minimum values as 2 σ uncertainty related to methodological choices. Our temperature model used in all optimizations is one-dimensional. This enables large numbers of computations and is appropriate because the surface elevation and temperature vary little upstream of the WDC site. The accumulation rate, however, increases upstream by ∼ 0.12 m ⋅ y − 1 in 30 km, and thus the vertical velocity field and temperature pattern are 2D. This variation is mostly accounted for automatically in our analysis because it is already embodied in the observed depth–age relation and hence in reconstructed accumulation rates, model velocities, and vertical heat advection. To examine the importance of residual 2D effects, we used a 2D flowline heat and ice flow model (62, 63) to calculate the depth–age and depth–temperature relations at the WDC site. Using these as inputs to the one-dimensional model and optimization, we compared temperature reconstructions for cases in which the upstream accumulation rate was either uniform or increasing as observed. LGM temperatures differed by as much as ± 0.5 °C for a variety of cases in which the imposed temperature history was not specified purely from the isotopic data. Such discrepancies represent another source of uncertainty on our reconstructed temperatures, and we treat this as a uniformly distributed random variable in the interval (−0.5, 0.5) °C, which contributes a 2 σ = ± 0.58 C ∘ to the limits shown in Fig. 2B (main text).

SI Influence of Elevation and Thickness Changes on Interpretation of Reconstructed Temperatures Our reconstruction is of temperature at the ice-sheet surface. Some aspects of the deglacial and Holocene temperature variations likely arise from changes of surface elevation via the atmospheric lapse rate, rather than from climate change connected to forcings and transports. This effect would warm the surface climate by about 0.7–1.0 °C per 100 m of subsidence. Further, the history of ice thickness influences our reconstructed surface temperature directly; thinning is associated with increased vertical ice velocity, which shifts reconstructed LGM temperatures to warmer values. Although the sensitivity depends importantly on the sequence and timing of thickness changes, a 130-m thinning from late glacial to mid-Holocene (or 100 m of surface subsidence after isostatic adjustment) raises the reconstructed LGM temperature by between roughly 0 °C and 0.3 °C, depending on the magnitude of preceding deglacial thickening. The lapse rate effect is thus the larger one. Following the LGM, the dominant glaciological process would have been thickening ice and a rising surface due to greatly increased accumulation rate, proceeding more rapidly at the higher-accumulation WAIS Divide than at the drier East Antarctic sites (64). In contrast, the rise of sea level between 18 kyBP and 8 kyBP imposed a relative elevation drop of about 140 m. And, during the Holocene, retreat of grounded ice in the Ross Sea (65) caused thinning and surface lowering. Geological evidence constraining this history is sparse. Cosmogenic exposure-age dates on nunataks at the margin of the ice sheet (64, 66) indicate the interior of the WAIS was thicker at 10 kyBP than present by < 100 m. The elevation histories from numerical full ice-sheet model simulations disagree by a few hundred meters despite use of similar climate histories. The model with the best treatment of ice dynamics of the grounded ice sheet (31, 59) yields thickness changes < 50 m from the LGM to present. Another recent model yields LGM thickness and elevation of ∼ 360 m and ∼ 260 m greater than present (60), although this has been calibrated to target an assumed thickness increase of 200 m. Two earlier models produced LGM thicknesses a few hundred meters thinner (30) and ∼ 100 m thicker (29) than present. All such simulations show significantly smaller elevation and thickness changes than suggested by the ICE-5G and ICE-6G glacial isostatic adjustment model scenarios commonly used in climate modeling experiments (67), which are almost certainly incompatible with both geological and ice-dynamical constraints. None of the ice-sheet models driven by climate records show more than 300 m of thinning at the WAIS Divide in the past 15 ky (29⇓–31). A change of this magnitude, if preceded by thickening due to the known accumulation increase, warms our reconstructed LGM temperature by ∼ 0.5 C ∘ . For this reason, the limits of our LGM temperatures shown in Fig. 2B (main text) include a term of this magnitude. Given the other factors contributing to uncertainty, this range expansion is also equivalent to a 2 σ temperature variation arising from 400-m thinning without prior thickening, an unlikely case. The preceding applies to the surface temperature reconstruction. For interpretations of constant-elevation climate, the larger lapse-rate effect must be included. As an outside case, suppose the ice at the WAIS Divide was 500 m thicker and the surface 400 m higher at LGM. Adjusting our constant-thickness LGM-to-late Holocene optimized temperature change of 11.3 ± 1.3 C ∘ with a + 1.5 C ∘ reconstruction effect and a +2.8–4 °C lapse rate effect would give a constant-elevation temperature change of only 6.4 ± 1.9 C ∘ . A more likely scenario of 200 m thickening when accumulation rises after the LGM, followed by 300 m of thinning to the mid-Holocene, would imply a + 0.53 C ∘ reconstruction effect and (assuming isostatic adjustment) a +0.5–0.7 °C lapse rate effect, giving a constant-elevation temperature change of 10.2 ± 1.4 C ∘ , which overlaps the range of surface temperatures shown in Fig. 2B. The difference between these cases illustrates that elevation changes could have a large influence on the interpretation of our temperature history. On the other hand, glaciological reasoning, ice sheet model simulations, and geological constraints all suggest that modest elevation changes are more likely, indicating that a smaller or possibly even negligible correction between surface temperature and constant-elevation temperature is more appropriate. In addition, the surface temperature difference between the LGM and the late Holocene is not much larger than that between the LGM and the early Holocene (Table S3), so elevation changes occurring after 10 kyBP have little effect on inferred deglacial warming. But if future evidence establishes that large ice-sheet changes occurred in the pre-Holocene, our interpretations related to constant-elevation climate will have to be revised.

SI Climate Forcings During Deglaciation At present the ocean transfers heat northward across the equator at a rate of about 1 PW (68). This amounts to an effective climate forcing throughout the Southern Hemisphere of approximately −4 W ⋅ m − 2 . In the reduced transport period (which culminated between 17 ka and 14.5 ka), the contribution to Antarctic warming would have been substantially less than 4 W ⋅ m − 2 , because the background circulation was weaker in the glacial climate, the transport did not cease entirely, and atmospheric heat exchange partly compensates (69⇓–71). Another contributor to West Antarctic warming was annually integrated insolation, which increased steadily through the deglacial period in the Southern Hemisphere. The increase between 18.5 ka and 15 ka, specifically, was ∼ 2.5 W ⋅ m − 2 at the periphery of Antarctica. The high albedo of the inland ice sheet ( ∼ 0.8) and moderate albedo of surrounding ocean regions ( ∼ 0.5) reduce the climate forcing to less than 1 W ⋅ m − 2 . The combined forcing from ocean transport and insolation changes was thus not much larger than 1 W ⋅ m − 2 . Other factors that might influence warming over the ice sheet include changes in the efficacy of atmospheric transport of heat and moisture and changes of winds at the Southern Ocean surface that alter sea surface temperatures around the continent. Such processes could be responding to variations of tropical climate or of ice-sheet topography in either hemisphere, for example. These hemisphere-specific forcings can be compared with the global forcings governing the difference between planetary temperatures of the LGM and Holocene. Taken together, decreasing ice-sheet extent (hence reduced albedo) and increasing greenhouse gas concentrations ( CO 2 and CH 4 ) account for about 90% of the direct global-averaged climate forcing, with respective magnitudes of ∼ 2.9 W ⋅ m − 2 and 2.2 W ⋅ m − 2 (2). (The very important variations of water vapor, snow cover, and sea ice constitute feedbacks rather than direct forcings.) Fig. 3B plots histories of these two forcings. To facilitate comparisons of relative timing, the histories are all scaled to range from zero to unity between 18.5 ka and 1 ka, times of nearly identical annually integrated insolation at 70 S ∘ . The ice-sheet albedo effect is concentrated in the Northern Hemisphere (72), and indeed this largest forcing appears to contribute little directly to the early deglacial changes in Antarctica. The Northern Hemisphere ice sheets probably instigated changes in the oceanic heat transport and sea ice through releases of meltwater and orographic effects on winds (73, 74), but the primary albedo forcing can be discounted; its magnitude by 15 ka was only ∼ 0.3 W ⋅ m − 2 . Antarctic temperature covaries much more strongly with the greenhouse gas forcing (41) (Fig. 3B, main text). Part of this correspondence likely arises because the climate warming itself and associated influences on Southern Ocean winds, sea ice, and currents drove the release of CO 2 and its atmospheric increase (75), and part arises because of the feedback from the increasing greenhouse effect. Antarctic warming leads the greenhouse gas forcing through the deglaciation, and the fractional increases (Fig. 3B, main text) put an upper limit on the contribution of the gases to Antarctic warming between 18.5 ka and 15 ka of ∼ 65%. The magnitude of the gas forcing by 15 ka was ∼ 1 W ⋅ m − 2 , compared with a forcing of 1 to a few W ⋅ m − 2 from combined insolation and reduced ocean transport. Thus, given this similarity, the observation that CO 2 could account for roughly half of the warming is consistent with normative estimates of its climatic impact.

Conclusion The temperature reconstruction reported here provides information about the thermal state of firn and ice through time at the West Antarctic Divide site, an essential input to studies of the depth–age relationship, interhemispheric climate phasing, CO 2 history, and covariation of temperature with accumulation (16, 47, 48). Of greatest immediate interest, however, is our demonstration that the global deglacial temperature change was amplified by a factor of 2–3 in the Antarctic, that Antarctic warming was largely achieved by 15 ka in coherence with records from Southern Hemisphere mountain ranges, and that climate models of the deglaciation perform well on average, but that the ones with lowest sensitivity can be discounted. The early warming of the Southern Hemisphere, which our study helps to quantify, arose from combined effects of reduced northward oceanic heat transport, increased insolation, and increasing atmospheric CO 2 . Quantitative simulation of this phenomenon could provide an illuminating challenge for model studies.

SI Measurement of Borehole Temperatures Temperatures in the fluid-filled WDC borehole were logged during December 2011 and again during December 2014, using the USGS PTLS. With the PTLS, a thermistor probe is lowered into the borehole at a constant rate while the sensor’s resistance and depth are continuously recorded. The resistance is later converted to temperature, using a special calibration function. Instrumental noise and temperature fluctuations due to borehole fluid convection are then removed using a wavelet analysis, and a deconvolution is performed to account for the finite response time of the probe. A full description of the measurement system, measurement uncertainties, and data processing techniques is available (50, 51). Uncertainties primarily arise from three sources: temperature-sensor calibration uncertainties and measurement-circuit drift, thermal fluctuations within the borehole due to convection of the borehole fluid, and the thermal disturbance caused by drilling processes. Except for the extreme ends of the borehole, the standard uncertainties due to these sources are less than 3.3 mK, 1.25 mK, and 4 mK, respectively. The quadrature sum provides a combined uncertainty of ∼ 5.3 mK for the temperature measurements in the fluid-filled portion (depths > 96 m) of the borehole. Above 96 m we used the temperature profile previously measured in a neighboring air-filled hole (52), shifted uniformly to match the values determined by the more accurate PTLS measurements.

SI Age–Depth Relation and Layer Thicknesses The age–depth relation and annual layer thicknesses were determined by identifying and counting annual layers back to ∼ 31 ka and by cross-correlations of gas records prior to this time (16, 47, 53). The former is completely independent of the reconstructed temperatures, whereas the latter uses them in the calculation of age offset between gas and ice. The 2 σ uncertainty on age offset at ages greater than 31 ka is less than ± 100 y. Optimizations using an age scale stretched or compressed by this amount change the reconstructed LGM temperature by a trivial amount, about ±0.05 °C. The range of our reconstructed temperature history (Fig. 2B, main text) includes the effects of stretching and compressing the age scale by 2 ky at 40 ka.

SI Comment on the Use of Nitrogen Isotopes Our use of the δ 15 N firn-thickness proxy in the reconstruction ( Δ T N in Eq. S3) is a unique application of a variable whose temperature dependence reflects a somewhat complex physical process. We emphasize that the optimization is still performed by use of the borehole temperatures alone. The improved model performance when Eq. S3 replaces Eq. S2 (Table S2) arises because the three main features of Δ T N (the mid-Holocene maximum, the warmer late glacial, and the colder LGM) are all favored independently by the borehole measurements. More precisely, if we exclude the most recent 2.5 ka when the prior model is already tightly constrained, then assigning a value ω = 1 in Eq. S3, meaning that Δ T N is included in accordance with theory, improves the borehole temperature match even though the optimization uses the same number of free parameters. This behavior provides evidence that the firn-thickness proxy serves as a thermometer despite potential complexities in the controls on firn structure and gas transport. The number of simultaneous free parameters in all optimizations (Eqs. S1–S3) remains constant (six).

SI Principal Features of Reconstructed West Antarctic Temperature Table S3 summarizes key quantitative features of our reconstruction.

SI Thickness History The forms of thickness histories H ( t ) used in optimizations are motivated by several considerations: (i) The doubling of accumulation rate between 18 ka and 15 ka would have caused thickening of a few hundred meters, (ii) retreat of the ice margins would have caused subsequent thinning, (iii) glacial geologic evidence indicates < 100 m of thinning since 10 ka, (iv) the correspondence of East and West Antarctic temperatures from 8 ka to the present indicates no significant elevation changes in this period, and (v) full ice-sheet model simulations yield LGM thicknesses relative to present ranging from about −200 m to +360 m. Our standard model scenario specifies 150 m of thickening between 18 ka and 15 ka followed by 50 m of thinning. In comparison, holding the thickness constant changes the reconstructed LGM temperature by +0.18 °C. Specifying a constant thickness before 12 ka, followed by thinning until 5 ka, changes the reconstructed LGM temperature by −1.59 °C to +1.42 °C as the initial thickness relative to modern ranges from −300 m to +450 m (Fig. S3). Including a 200-m increase between 18 ka and 13 ka, before the 12-ka to 5-ka thinning, reduces the LGM temperature by ∼ 0.5 °C compared with the previous case. Using thickness histories calculated for the WDC site in specific ice-sheet model simulations yields the following changes to LGM temperature: −0.085 °C [model of Pollard et al. (59)], +0.55 °C [model of Golledge et al. (60)], and +0.74 °C [model of Pollard and DeConto (29)].

Acknowledgments We are deeply indebted to many participants in the WDC project and especially thank K. Taylor, E. J. Brook, and J. W. C. White. The helpful comments of two anonymous referees are gratefully acknowledged. This work is funded through the US National Science Foundation Grants 0539232, 0537661 (to K.M.C.), 0537930, 1043092 (to E.J.S.), 1043518 (to C.B.), 0944199, 0944197, 0440666 (to E.D.W.), 0539578, 1043528, 1338832 (to R.B.A.), and 0538657 (to J.P.S.) and National Aeronautics and Space Administration Grant NNX12AB74G (to M.K.). We gratefully acknowledge additional support from the Martin Family Foundation (K.M.C.), the USGS Climate and Land Use Change Program (G.D.C.), and National Oceanic and Atmospheric Administration Climate and Global Change Fellowships (C.B.).

Footnotes Author contributions: K.M.C. and C.B. designed research; K.M.C., G.D.C., E.J.S., C.B., T.J.F., M.K., E.D.W., R.B.A., and J.P.S. performed research; K.M.C., G.D.C., E.J.S., C.B., T.J.F., M.K., and E.D.W. analyzed data; and K.M.C. wrote the paper.

The authors declare no conflict of interest.

This article is a PNAS Direct Submission.

This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10.1073/pnas.1609132113/-/DCSupplemental.