Articles | Volume 18, issue 6
Research article
21 Jun 2022
Research article |  | 21 Jun 2022

Improving temperature reconstructions from ice-core water-isotope records

Bradley R. Markle and Eric J. Steig

Oxygen and hydrogen isotope ratios in polar precipitation are widely used as proxies for local temperature. In combination, oxygen and hydrogen isotope ratios also provide information on sea surface temperature at the oceanic moisture source locations where polar precipitation originates. Temperature reconstructions obtained from ice-core records generally rely on linear approximations of the relationships among local temperature, source temperature, and water-isotope values. However, there are important nonlinearities that significantly affect such reconstructions, particularly for source region temperatures. Here, we describe a relatively simple water-isotope distillation model and a novel temperature reconstruction method that accounts for these nonlinearities. Further, we examine in detail many of the parameters, assumptions, and uncertainties that underlie water-isotope distillation models and their influence on these temperature reconstructions. We provide new reconstructions of absolute surface temperature, condensation temperature, and source region evaporation temperature for all long Antarctic ice-core records for which the necessary data are available. These reconstructions differ from previous estimates due to both our new model and reconstruction technique, the influence of which is investigated directly. We also provide thorough uncertainty estimates for all temperature histories. Our reconstructions constrain the pattern and magnitude of polar amplification in the past and reveal asymmetries in the temperature histories of East and West Antarctica.

1 Introduction

Stable-isotope ratios of water have been the foundational proxy for polar paleoclimate research for more than a half-century (Langway1958; Gonfiantini1959; Dansgaard1964). Primarily used as a temperature proxy, stratigraphic records of water-isotope ratios in ice sheets provide detailed histories of Earth's climate over hundreds of thousands of years (Dansgaard et al.1969; Petit et al.1999), providing insight into the past magnitudes, spatial patterns, and phasing of climate change across the globe (Masson-Delmotte et al.2006; EPICA Community Members2006; WAIS Divide Project Members2013, 2015). Both oxygen and hydrogen have stable isotopes whose ratios (18O/16O and 2H/1H) are commonly expressed as deviations, δ18O and δD, from Vienna Standard Mean Ocean Water (VSMOW) in per mill (‰):

(1) δ = R x - R std R std ,

where Rx is the ratio in the sample and Rstd is the ratio in VSMOW.

Poleward transport of moisture by the climate system, the progressive removal of moisture from the atmosphere by condensation and precipitation, and the fractionation of water-isotope ratios during phase changes are all processes inherently linked to temperature and together underpin the use of water-isotope ratios in polar precipitation as a temperature proxy (Craig1961; Epstein et al.1963; Dansgaard1964; Gonfiantini1965). The strong empirical correlation between the water-isotope ratios in precipitation and surface temperature supports this interpretation (Petit et al.1999; Jouzel et al.1997; Masson-Delmotte et al.2008). Air temperatures during condensation (Petit et al.1999; Jouzel et al.1997) and during initial moisture evaporation (Vimeux et al.2002) can be reconstructed from ice-core water-isotope records if the relevant scaling relationships can be determined from theory, models, or observations (Vimeux et al.2002; Kavanaugh and Cuffey2002; Stenni et al.2010). Here, we examine the widely used assumption of linearity in the scaling relationships between water-isotope ratios and temperature.

1.1 Temperature reconstructions

Any interpretation of water-isotope ratios as a proxy for temperature requires a model, whether conceptual, statistical, or numerical. A conceptual model of progressive distillation and integrated fractionation (e.g., Dansgaard1964) is sufficient to qualitatively interpret variations in water-isotope ratios as variations in temperature in the high latitudes. The simplest quantitative interpretation of ice-core water-isotope records relies on the empirical correlation between observed water-isotope ratios of precipitation and surface temperature at the precipitation site (Petit et al.1999; Jouzel et al.1997). A limitation of this approach is the possibility to conflate the “spatial slope” between water isotopes and temperature, which is the relationship observed across a range of modern sites, and the “temporal slope”, which is the relationship at a single point through time (Jouzel et al.1997). Relatedly, this approach also does not account for simultaneous and independent changes in evaporation conditions, which can impact high-latitude water-isotope ratios in several ways. Initial evaporation temperature, together with the condensation temperature, determines the total temperature gradient through which moisture must be distilled to reach a given site. Further, evaporative conditions set the initial isotopic values of the vapor before distillation. The isotope ratios of vapor above the ocean depend on the temperature during evaporation, the isotopic values of the seawater, and the occurrence of kinetic fractionation during evaporation, which is driven by sub-equilibrium relative humidity and influenced by sea surface temperature and wind speed (Merlivat and Jouzel1979; Jouzel et al.1982).

A more complete approach to reconstructing temperature from water-isotope records is to employ numerical models that account for the combined influence of variability in both evaporation and condensation temperatures, as well as other factors. Reconstructing two unknowns (i.e., both evaporation source and condensation site temperatures) requires two constraints, which are provided by the oxygen and hydrogen isotope ratios and the relationship between them. While the oxygen and hydrogen isotope systems have similar behavior in the atmosphere, there are differences in their response to the same environmental conditions and to processes such as kinetic fractionation. The deuterium excess is the weighted difference between δ18O and δD, dxs=δD-8×δ18O, and is commonly used to quantify these differences (Dansgaard1964; Merlivat and Jouzel1979).

Changes in water-isotope parameters measured in precipitation at an ice-core site, Δδ18O and Δdxs, can be conceptualized as driven by changes in site and evaporation source temperature, ΔTsite and ΔTsource:


where β and γ are the partial derivatives of δ18O and dxs with site and source temperature, respectively. The magnitudes of β and γ can be diagnosed from water-isotope distillation models for the ice-core site in question (Vimeux et al.2002; Kavanaugh and Cuffey2002; Stenni et al.2010; Uemura et al.2012). Once these slopes are established, the equations may be solved for ΔTsite and ΔTsource using records of δ18O and dxs (Vimeux et al.2002; Stenni et al.2010; Uemura et al.2012).

1.2 Nonlinearities in isotope fractionation and the deuterium excess definition

The temperature reconstruction approach described above depends on the assumption that the parameters, β and γ, are fixed in time and independent of temperature. However, the β and γ parameters, as diagnosed from model simulations, are found to be different for different ice-core sites, whose only distinguishing characteristics are differing modern surface conditions (e.g., Stenni et al.2010; Uemura et al.2012). Thus, β and γ depend on the site conditions, which obviously change over time.

Another issue with the linear reconstruction approach is the definition of the deuterium excess parameter (Uemura et al.2012; Markle et al.2017). The origin of the slope in the definition of deuterium excess is an empirical fit to global precipitation measurements (Dansgaard1964). However, a linear relationship between δ18O and δD is not fundamental (Craig1961); equilibrium fractionation alone drives a nonlinear relationship between δ18O and δD (Markle et al.2017). While the effects of source region conditions on deuterium excess of vapor are nearly linear during initial evaporation (Merlivat and Jouzel1979; Uemura et al.2008), the signal is not uniformly preserved as moisture is transported toward the deposition site. Kinetic fractionation that occurs during transport (Jouzel et al.1982) alters the deuterium excess of the vapor, as does equilibrium fractionation during condensation, owing to biases in the linear definition (Markle et al.2017). Thus, the sensitivity of dxs in precipitation to evaporation and condensation temperatures must vary as a function of the total condensation and fractionation experienced during transport to any deposition site and is thus a function of Tsite.

Some of these issues have been addressed by redefining the deuterium excess parameter (Uemura et al.2012; Markle et al.2017). Uemura et al. (2012) fit a second-order polynomial to a compilation of δ18O and δD data, where δx=ln1+δx, and defined a phenomenological, nonlinear deuterium excess parameter:

(4) d ln = δ D - A × ( δ 18 O ) 2 + B × δ 18 O ,

with coefficients A=-28.5 and B=8.47 (note that the coefficients and δ values are unitless; for example, δ18O=0.040 not 40, with the ‰ ignored).

While still empirical, this definition of deuterium excess reduces the influence of kinetic fractionation during transport and the biases inherent to the linear definition, making it a more faithful qualitative proxy for source region conditions (Uemura et al.2012; Markle et al.2017), and is particularly important at the coldest Antarctic sites where nonlinear effects overwhelm the dxs definition. However, the same distillation processes that lead to biases in the linear definition of the deuterium excess parameter will also bias the results of temperature reconstructions if fixed sensitivities (Eqs. 2 and 3) are assumed.

Here we examine these issues in water-isotope-based temperature reconstructions and suggest an improved technique.

2 A (relatively) simple water-isotope model

The quantitative reconstruction of temperatures from water-isotope ratios rests on the encapsulation of fractionation processes in models. Any investigation into nonlinearity in those relationships will depend on the representation of those physics. To assess the importance of those nonlinearities, we construct a model that is relatively simple while still faithfully representing the observed relationships between the hydrogen and oxygen isotope ratios in polar precipitation. We describe the construction of the Simple Water Isotope Model (SWIM) in detail in the Appendix. Here we describe the conceptual framework of the model.

The underpinning of SWIM is shared by many water-isotope models: the transport and distillation of moisture down climatological temperature gradients. Moisture is evaporated from the oceans in the low and midlatitudes and transported toward the poles. As air cools, the saturated vapor pressure decreases nonlinearly, and moisture above saturation is removed by precipitation. During these phase changes, water fractionates; the vapor and precipitation falling from it become increasingly depleted in the heavier isotope. The total fractionation at any point is a consequence of the temperature gradient through which the water is distilled, as well as the mean temperature of that gradient, owing to nonlinearity in the Clausius–Clapeyron relationship. A change in the average condensation temperature at a site thus results in a change in the isotope ratios of precipitation at that site. This is the essential (though not sole) reason that high-latitude water-isotope ratios are a useful temperature proxy. They are driven by two basic nonlinear processes, the Clausius–Clapeyron relationship and Rayleigh distillation (see Sect. A2.2).

Other processes can be important as well. The temperature dependence of fractionation factors, for example, generally amplifies the temperature relationship. While any single precipitation event at a site may be subject to a variety of additional factors and processes, the long-term mean is strongly influenced by climatological moisture distillation.1

Our model distills moisture down thermodynamic pathways defined by temperature and pressure. Changes in water-isotope ratios are driven neither by changes in space nor time but by changes in the thermodynamic variables that cause the water to change phase. The temperature gradient of the pathway is prescribed from an initial evaporation temperature, T0, to a final condensation temperature, Tc. The pathways are pseudo-adiabatic, consistent with isentropic moisture transport to the Antarctic (Bailey et al.2019) and the basic assumption of Raleigh distillation that moisture is removed after precipitation. A superposition of many thermodynamic pathways is required to represent a single Antarctic precipitation site, reflecting both the range of precipitation conditions experienced at a site and moisture transport from sources with a distribution of evaporative conditions (Markle et al.2017, Fig. 1). An example of a set of these pathways is shown in Fig. 2. We use climatological correlations to relate initial evaporation air temperature, T0, to other initial conditions including sea surface temperature, SST0, and relative humidity, RH0 (see Sect. A1.1).

Figure 1Moisture sources and transport to Antarctica from the moisture-tagged Community Atmosphere Model (CAM) experiment (Markle et al.2017). (a) Difference, in degrees latitude, between the latitude of precipitation and mean latitude of evaporation (effectively mean transport in degrees of latitude to any site). (b) Mean latitude of evaporation vs. latitude of precipitation. All longitudes are shown by the black line; longitudes encompassing West Antarctica are shown by the blue line, while longitudes encompassing East Antarctica are shown in red. (c) The mean evaporation latitude (in  S) of precipitation falling at all Antarctic grid points. Select ice-core sites shown in white. (d) The relationship between mean evaporation latitude and site elevation across Antarctica. Select ice-core locations shown in color; see text for site details. (e) Latitudinal moisture source distributions for select Antarctic ice-core sites, colored by site elevation.

Figure 2Simple Water Isotope Model. (a) Condensation temperature–pressure pathways for pseudo-adiabatic distillation. Pathways are colored by initial evaporation air temperature in all panels. (b) Mixing ratios for all pathways. (c) δ18O of precipitation for all pathways. (d) Relationship between δ18O and dln of precipitation for all pathways. Black dots show water-isotope values of annual averaged precipitation from a compilation of global observations (see text), while grey dots show a broader set of monthly observations from the GNIP database (IAEA2001).


We consider the temperature dependence of kinetic and equilibrium fractionation during both evaporation and transport, as well as mixed liquid and ice phases of precipitation. The model incorporates supersaturation at very cold conditions, which is tuned to match the observed relationship between the oxygen and hydrogen isotope ratios in global precipitation (Fig. 2d) rather than the relationships between those parameters and climate variables such as temperature. We investigate the sensitivity of the model and the resulting reconstructions to uncertainties and assumptions including fractionation factors, evaporative closure assumptions, precipitation schemes, supersaturation, the pseudo-adiabatic pathway, initial climatological correlations during evaporation, and non-uniqueness. We also investigate the influence of mixing during both evaporation and transport, the potential influence of seasonality (and intermittency) of precipitation and evaporation, and the relationship between surface temperature at a site and the vertically integrated, precipitation-weighted condensation temperature above that site. We make no attempt to model post-depositional processes.

3 Temperature-dependent slopes

To investigate the sensitivity of water-isotope ratios of Antarctic precipitation to site and source conditions, we use SWIM to model isotopic state spaces. We run SWIM through a large ensemble of temperature pathways defined by T=T0,T0-dT,,Tc, with dT=0.1C. We run nearly 24 000 trajectories to fill the plausible parameter space of 0 C T028C and −70C Tc10C. We first examine the model with base assumptions and parameterizations2. We next investigate the sensitivity to these choices.

The modeled isotopic state spaces are shown as maps whose x and y coordinates are the condensation and evaporation temperatures, Tc and T0, respectively, and whose z (color) dimension is the isotopic value of final precipitation in Antarctica: δ18O in Fig. 3a, δD in Fig. 3b, dxs in Fig. 3c, and dln in Fig. 3d.

Figure 3Water-isotope state spaces as a function of the boundary conditions Tc, the condensation temperature, and T0, the mean evaporation temperature. Surface shading and contour lines are the water-isotope values of precipitation (in units of ‰) (a) δ18O, (b) δD, (c) linear deuterium excess dxs, and (d) nonlinear deuterium excess dln.


The gradients of both the δ18O and δD surfaces are predominantly in the direction of the condensation temperature (the x axis in Fig. 3), emphasizing the strong condensation temperature dependence of these parameters. However, the slopes of both δ18O and δD are not strictly linear with condensation temperature Tc, clearly varying with its absolute value (and to a much lesser extent with the evaporation temperature, T0, due to its influence on the total distillation gradient). Further, the partial slopes of δ18O and δD with respect to the evaporation source temperature depend strongly on the absolute values of both the evaporation and condensation temperatures, evidenced by the changing angle of the contour lines in Fig. 3. The partial derivatives of the isotopic surfaces with respect to T0 and Tc are shown explicitly in Figs. 4 and 5. It is important to recognize that the partial derivatives with respect to T0 are not for the initial vapor at the point of evaporation, but for the precipitation after the vapor has passed through the distillation pathway to the final precipitation site. The sensitivity of isotopic values of precipitation to source region conditions is a function of the total distillation that the moisture experiences.

Figure 4Partial derivatives of isotope state spaces with respect to condensation temperature: (a) δ18OTc, (b) δDTc, (c) dxsTc, (d) dlnTc. Shading and contours in all panels represent the slope in ‰ C−1.


Figure 5Partial derivatives of isotope state spaces with respect to evaporation temperature: (a) δ18OT0, (b) δDT0, (c) dxsT0, (d) dlnT0. Shading and contours in all panels represent the slope in ‰ C−1.


The modeled dxs surface shows strong slopes along both the condensation temperature and evaporation temperature axes (Fig. 3c), as does modeled dln (Fig. 3d). The dln depends more strongly on the evaporation temperature than the dxs. In particular, at the coldest condensation temperatures, variability in dxs is dominated by the condensation temperature, reflecting the influence of kinetic fractionation during condensation and the nonlinear bias inherent to the historical linear definition (Uemura et al.2012; Markle et al.2017). These model results (Figs. 3, 4, and 5) demonstrate that the logarithmic definition of the deuterium excess parameter (dln) is a more faithful qualitative proxy for source region conditions than the linear definition (dxs). Even at very low condensation temperatures, dln still depends strongly on the initial evaporation temperature, whereas linear dxs becomes more dependent on condensation temperature.

A “bump” in the partial derivatives of all isotope parameters with respect to the condensation temperature is seen around −30C, arising primarily from the transition between liquid and ice condensate (Fig. 4), whose relationship to temperature is prescribed in the model and based on satellite data (see Appendix A2.2). The slopes in this region also depend on the parameterization of supersaturation. This local change in the partial derivatives with respect to Tc is smoothed somewhat when atmospheric mixing is incorporated into the model (see Appendix A5). The changes in the partial derivatives with respect to Tc across the entire range investigated are larger than these changes localized around −30C.

The temperature reconstruction technique described in Sect. 1.1 is based on the linearization of the slopes between δ18O, dxs, condensation temperature, and evaporation temperature, and it assumes that the β and γ parameters in Eqs. (2) and (3) are fixed over the range of reconstructed ΔTsite and ΔTsource. Our results demonstrate that the assumption that β and γ are independent of temperature (i.e., that the sensitivities are fixed and linear) is problematic. The parameters γ1 and γ2 in Eq. (2) are comparable to the slopes δ18OTc and δ18OT0 in Figs. 4 and 5. Although the slope of δ18O along the condensation temperature axis, δ18OTc, does not change dramatically, it is clearly variable, as is δ18OT0. The slopes dxsTc and dxsT0 (comparable to β1 and β2 in Eq. 3, respectively) are highly variable. Indeed, the dxs surface in Fig. 3 has a saddle at moderate condensation temperatures, over which dxsTc changes sign (Fig. 4c). This shows that the assumption of constant β and γ parameters in isotope-based temperature reconstructions is incorrect and may be reasonable only under narrow ranges of ΔTsite and ΔTsource. For plausible changes in site temperatures, assuming a fixed γ, for example, may not only lead to errors in magnitude but even to errors in the sign of γ and ultimately ΔT. The use of a fixed γ can introduce spurious variability into temperature reconstructions, particularly of T0, the evaporation source temperature.

Evidence for the critical change in the sign of dxsTc (Fig. 4c) can actually be seen directly in Antarctic ice-core records. Deuterium excess (dxs) records from core sites whose average conditions lie on the same side of this change in slope are generally correlated over the last 60 000 years, while sites whose average conditions lie on opposite sides of the change in slope are weakly correlated or even anticorrelated with each other (Fig. 6). Note how glacial–interglacial changes in dxs at the coldest sites (blue colors in Fig. 6c) are opposite to those at the warmest sites (red and orange colors). On the other hand, dln values from different cores are vastly more coherent (Fig. 6e and f) for reasons described above. Comparison of dxs and dln between cores (Figs. 6 and A24) together with our modeling results demonstrates that the incoherence of Antarctic dxs records arises from the change in slope of dxsTc, not from differing source conditions between sites. This exposes a fundamental flaw in the assumptions used in traditional isotope temperature reconstructions.

Figure 6Time series and cross-correlation matrices for eight different deep Antarctic ice-core sites. (a, b) δ18O; (c, d) dxs; and (e, f) dln. In the time series plots, each record is colored by its Holocene average δ18O value. See Sect. 4.2 in the text for details and references for the ice-core records. All original records are interpolated to even 50-year time spacing on the Buizert et al. (2018) synchronized timescale where possible, or they are plotted on original published timescales. All records are ordered by their approximate modern surface temperature in the cross-correlation matrices.


The issue of variable isotope–temperature scalings is implicit in previous work. Uemura et al. (2012), for example, following Stenni et al. (2010), used an isotope model to calculate the relevant β and γ parameters for several East Antarctic ice-core sites. Using the same isotope model, they calculated different scalings for each site. However, by assuming that these slopes are constant for each site, they do not consider the possibility that one site's conditions may have been more like another's in the past. Recognizing this as well as the inability of their model to simultaneously match observed site temperature and δ values, Uemura et al. (2012) create several reconstructions for the Dome Fuji site utilizing different linearizations of the model. They do not, however, attempt a reconstruction that accounts for the nonlinearities in the water-isotope–temperature relationships.

The solution to this issue of slope nonlinearity, within the linear isotope temperature reconstruction framework (Eqs. 2 and 3), is not obvious. The nonlinearities in the slopes of the isotope surfaces depend on the absolute condensation and evaporation temperatures that represent the target of the reconstruction, which are of course not known a priori. We next present a novel temperature reconstruction framework, which takes into account the inherent nonlinearities in the water-isotope fractionation process.

4 Nonlinear temperature reconstructions

4.1 Reconstruction method

For every pair of T0 and Tc inputs to SWIM there is a corresponding modeled value of δ18O, δD, and dln of final precipitation as shown in Fig. 3. We invert the modeled state spaces and project each independent temperature parameter onto a pair of dependent isotope values, e.g., δ18O and dln. This defines a set of maps, with x and y axes of δ18O and dln and with z axes Tc and T0, as shown in Fig. 7. To reconstruct Tc and T0, the inverted model results may be used as a lookup table: a pair of δ18O and dln measurements determine a pair of Tc and T0 reconstructions (Fig. A16). While previous reconstruction methods (e.g., Vimeux et al.2002; Kavanaugh and Cuffey2002; Stenni et al.2010) linearize the slopes calculated by a water-isotope model around the modern climate state, this method accounts for the changes in the slopes that depend on the mean state. Further, there is no need to find analytical solutions to the model or fit families of high-order polynomials to the results.

Figure 7Inverted T0 and Tc surfaces as a function of modeled dln and δ18O of precipitation. (a) Surface shading and contours show condensation site temperature, Tc, in degrees Celsius (C) as a function of dln and δ18O of precipitation at a site. (b) Surface shading and contours show evaporation source temperature, T0, in degrees Celsius (C) as a function of dln and δ18O of precipitation at a site.


The boundary conditions Tc and T0 may be projected onto axes defined by any two isotope parameters, which may then be used to reconstruct temperature. Since the only unique isotope information comes from the original δ18O, δD measurements (dxs and dln being second-order parameters), any combination of those parameters may seem equally well suited for the purposes of temperature reconstruction. In practice, however, δ18O and dln represent the optimal pair of parameters to use for temperature reconstruction. This result is examined in more detail in the Appendix (see Sect. A6 and Fig. A15). The fundamental reason is that the logarithmic excess parameter, as a second constraint, provides an axis more orthogonal to the variability we are attempting to reconstruct. This is also the same reason dln is a better qualitative proxy for source region temperature than dxs. After proposing the dln parameter, Uemura et al. (2012) suggest that there is no added value in the logarithmic parameter over the traditional linear dxs in terms of the temperature reconstruction equations. While true for the linear temperature reconstruction equations, this is not the case when the nonlinearities of β and γ are accounted for.

Note that we could just as readily reconstruct source region relative humidity (RH0) instead of source region air temperature (T0), since we assume climatological relationships between them. Indeed, RH0 is a strong lever on the deuterium excess of evaporation. We reconstruct T0 out of interest in the parameter from a climate dynamics perspective. In principle, our method can be extended to independently reconstruct more than two variables simultaneously, such as Tc, T0, and RH0, by modeling multidimensional parameter spaces. This of course requires additional measured constraints, which could include the δ17O of precipitation, the accumulation rate, or other variables modeled in this framework, and is discussed in the Appendix.

4.2 Absolute temperature reconstructions

An advantage of the reconstruction technique presented here is that we are able to reconstruct absolute evaporation and condensation temperatures, not just relative variability, as in previous techniques. There are several additional considerations that are important in making these reconstructions. First, we are often interested in the surface air temperature for paleoclimate studies rather than the condensation temperature. In Appendix A3.2 we review previous work on the relationship between surface and condensation temperature over Antarctica and conduct novel analysis using the high-resolution MERRA-2 reanalysis data (Gelaro et al.2017). We also examine the seasonality and vertical distribution of Antarctic precipitation and reevaporation. Based on these analyses we use a simple, linear temperature-dependent relationship to estimate weighted, annual mean surface air temperature from our condensation temperature reconstructions and account for the uncertainty in this relationship.

We examine in detail how assumptions and modeling choices concerning initial evaporation (Appendix A2.1) impact our results, as well as the potential for bias arising from seasonality in evaporation from the ocean. We also examine the influence of non-uniqueness on our temperature reconstruction technique arising from below-freezing evaporation conditions (Appendix A9.4).

Finally, we conduct an extensive uncertainty analysis of our temperature reconstructions (Appendix A9). We calculate numerous isotope state spaces from the same temperature parameter space using multiple iterations of the model in which model assumptions are altered and parameters are varied over plausible ranges. We calculate the uncertainty in our reconstructions arising from the supersaturation parameterization, the evaporation fractionation factors, the evaporation closure assumption, the precipitation scheme, the influence of vapor mixing during transport, and other model choices (see, e.g., Fig. A25). We use the ensemble of isotope state spaces to estimate both the absolute and relative uncertainty in our temperature reconstructions. As an example, the central reconstructions and their uncertainties for the evaporation and condensation temperatures for the WAIS Divide ice core (WDC) are shown in Fig. 8. Our reconstruction of relative temperature variability has much lower uncertainty than the reconstruction of absolute temperature, and we find lower uncertainty in the reconstruction of condensation temperature than evaporation temperature.

Figure 8Temperature reconstructions and uncertainty estimates for the WAIS Divide ice-core (WDC) site. (a) moisture source mean evaporation temperature and (b) ice-core site condensation temperature. Dark blue lines are the central estimate, while red lines show the bounds of relative temperature uncertainty and cyan lines show the bounds of the absolute temperature uncertainty. Reconstructions resampled to even 40-year spacing.


Figure 9Temperature reconstructions for seven ice-core sites: WDC, EDC, EDML, Siple Dome, Vostok, Dome Fuji, Talos Dome, and the South Pole ice core (SP). (a) Ice-core site locations on the Antarctic continent. Reconstructions of (b) moisture source evaporation air temperature, (c) ice-core site condensation air temperature, and (d) ice-core site surface air temperature. All records are resampled to even 200-year resolution for visual clarity. Light shading in panels (b–d) is the absolute temperature uncertainty, while dark shading shows the relative temperature uncertainty, though the latter can be difficult to see at the scale plotted here given the larger magnitude of point-to-point variability. For more detail see Figs. A17A18, and  A19.

We reconstruct condensation site and surface temperatures and evaporation source temperatures for eight different Antarctic deep ice-core sites for which there are δ18O and dln records (Fig. 9). The records include WDC (Markle et al.2017; WAIS Divide Project Members2013; Steig et al.2013) and Siple Dome (Brook et al.2005; Schilla2007) from West Antarctica, as well as the EDML (Stenni et al.2010), EDC (Stenni et al.2010), Vostok (Vimeux et al.2002), Dome Fuji (Uemura et al.2012), Talos Dome (Stenni et al.2011), and South Pole (SP, Steig et al.2021) records from East Antarctica. We correct all records for changes in the isotopic composition of seawater (Bintanja and Van de Wal2008), δ18Osw, following the method outlined in Uemura et al. (2012) and Stenni et al. (2010).

We make no corrections to the records for elevation changes or ice flow, lacking sufficient constraints for all records. These effects are likely small in both East Antarctica (Stenni et al.2011) and West Antarctica (WAIS Divide Project Members2013). Consideration of ice sheet changes are of course critical to disentangling the sources of temperature variability across timescales (e.g., Werner et al.2018) and require additional information and assumptions. However, for logical consistency our aim here is to reconstruct temperature as it has influenced water-isotope ratios recorded in the ice cores, leaving the identification of the causes of those changes in temperature for future work. This means that neither the sites of precipitation nor evaporation are strictly fixed in space. Indeed, variability in moisture source locations may cause a significant amount of variation in reconstructed evaporation temperature (Markle et al.2017). We aim to maintain the broad utility of our reconstructions by building in as few assumptions about the sources of the temperature variability as possible.

4.3 Comparison of linear and nonlinear reconstruction techniques

Linear temperature reconstruction using SWIM

We evaluate the significance of our approach by comparing our nonlinear reconstructions of condensation and evaporation temperature with reconstructions following the traditional linear approach. We calculate the equivalent linear β and γ coefficients for each of eight ice-core sites by regressing the Tc and T0 temperature fields to subsets of δ18O and dxs SWIM results representative of the Holocene (<10 ka) and Last Glacial Period (20 to 30 ka) intervals of each core. We find that the β and γ values significantly differ for each core and between the different time intervals owing to the temperature dependence of the sensitivities. In particular, β2 changes substantially between the glacial and the Holocene.

We reconstruct relative changes in moisture source and site surface temperature for each ice-core location using the two sets of linear β and γ coefficients found above and the traditional linear method. We then compare these linear reconstructions to our full nonlinear reconstruction and show the difference in reconstructed surface temperature and evaporation temperature in Fig. 10. We also show the residuals between the two linear reconstructions. In general, the residuals are not constant offsets but vary as a function of reconstructed temperature, demonstrating the temperature dependence of the slopes. They also differ, sometimes in sign, between the ice-core sites (Fig. A20). Linearization can obscure true variability or introduce spurious variability into the reconstructions, depending on the actual conditions of the site over time.

In the case of the surface temperature reconstruction, the errors introduced by linearization can be up to ±1C depending on the core site and are generally smaller for the colder East Antarctic sites. In the case of evaporation temperatures, the introduced errors are considerably larger with values up to ±2C. Further, the total variability in reconstructed evaporation temperature is much smaller than that in ice-core site surface temperature. The errors introduced into the reconstructed evaporation temperatures by ignoring the nonlinearities can be nearly as large as the total reconstructed variability. It is thus problematic to attempt reconstructing evaporation temperatures without accounting for nonlinearities.

In Appendix A8 we show that using the nonlinear reconstruction technique yields greater correlation amongst all records of Ts and especially of T0 (with increases in shared variance up to 38 %, Fig. A22). These results suggests that linear reconstructions have obscured coherent underlying climate signals, especially in evaporation temperatures. This same reasoning supports the qualitative use of the dln parameter over the linear dxs parameter (Fig. 6; Markle et al.2017).

Figure 10Differences between reconstructed (a) Ts and (b) T0 using different reconstruction techniques for multiple core sites. Blue lines show the difference between our full nonlinear reconstruction and a linear reconstruction using β and γ slopes linearized around Holocene conditions. Red lines show the difference between our full nonlinear reconstruction and a linear reconstruction using slopes linearized around glacial conditions. Purple lines show the difference between the two linear reconstructions.


The relationship between δ18O and Tc is largely linear across a wide range of values of Tc, regardless of evaporation temperature. The ice-core site temperature reconstructions from the linear and nonlinear reconstruction techniques have relatively small differences. However, as seen in Fig. 10, there are small artifacts arising from slight nonlinearity in the δ18O-to-Tc relationship, particularly for relatively warm sites like those in West Antarctica. The primary source of this nonlinearity is the change in total fractionation factor as the air parcel transitions between liquid-only and ice-only condensate. SWIM retains liquid condensate at colder temperatures than previous models (e.g., Kavanaugh and Cuffey2002), in line with satellite measurements (Hu et al.2010). The resulting transition of fractionation factors drives the nonlinearities in the δ18O-to-Tc relationship at temperatures relevant to West Antarctica, ultimately resulting in larger differences between the linear and nonlinear reconstruction techniques at those sites compared to cores from East Antarctica. Because our model uses a consistent supersaturation parameterization in the model's isotope and precipitation schemes, the relationship between δ18O and Tc is actually more linear in SWIM than in other comparable models.

In Appendix A10 we compare our temperature reconstructions for several East Antarctic ice cores to previously published reconstructions using the linear technique with coefficients estimated from different water-isotope models (Stenni et al.2004, 2010; Uemura et al.2012). The differences between our reconstructions and previous reconstructions arise from differences in both the reconstruction technique and the underlying isotope models. In general, the previously published linear reconstructions overestimate changes in both site and source temperature compared to our nonlinear reconstructions (Fig. 11). For example, previous reconstructions find larger warming since the Last Glacial Period for East Antarctic sites compared to ours (up to 18 % larger in the case of site temperature and up to 125 % larger for source temperatures). However, the differences between our reconstructions and previous reconstructions vary between sites owing to both temperature dependence and model differences, and in some cases they are quite small, such as with the Stenni et al. (2010) reconstructions. The residuals between relative temperature change in the nonlinear and previous linear reconstructions are shown in Fig. 12. The residuals are not random but rather correlated with the reconstructions themselves, reflecting the nonlinear biases. In general we find the largest differences, as a fraction of total reconstructed variability, in the moisture source temperature reconstructions. Indeed, these residuals can be over twice as large as the total reconstructed variability in moisture source temperature.

Figure 11Reconstructions of relative change in Antarctic surface temperature (ΔTsite, left panels) and source region evaporation temperature (ΔTsource, right panels) for four East Antarctic ice-core site: Vostok, EDC, EDML, and Dome Fuji. The nonlinear reconstructions (this study) are shown in black, while published linear reconstructions are shown for each site in color. The linear coefficients for the published reconstructions are compiled in Uemura et al. (2012) (see Tables 1 and 2). Linear methods labeled U12 for Vostok, EDC, and EDML were calculated by a simple Rayleigh-type model (Uemura et al.2012). Reconstructions U12a–e for Dome Fuji represent a sensitivity study from Uemura et al. (2012). Reconstructions S03 and S09 are from Stenni et al. (2004, 2010).

Figure 12Differences between our full nonlinear reconstruction and multiple previously published linear reconstructions (Stenni et al.2004, 2010; Uemura et al.2012) of (a) ice-core site surface temperature, (b) site condensation temperature, and (c) evaporation source temperature for multiple core sites. Full reconstructions are shown in Fig. 11.


5 Discussion

Using the self-consistent, nonlinear temperature reconstruction technique for eight different ice-core sites, we next investigate the patterns of Southern Hemisphere temperature change through time. In Fig. 9 we show the nonlinear reconstructions of Antarctic surface temperature and moisture source evaporation temperature for the eight ice-core records. At the WDC site in West Antarctica there is an independent estimate for the magnitude of glacial–interglacial temperature change from the borehole temperature profile (Cuffey et al.2016). Our results are in good agreement with those findings in both the absolute value of reconstructed temperatures and the magnitude of glacial–interglacial change. Cuffey et al. (2016) find 11.3±1.8C warming at WDC during the deglaciation; we reconstruct 11.2±0.5C of warming (calculated as the difference between the average surface temperature at 27–24 and 11–9 ka, for direct comparison to Cuffey et al.2016). Using an independent temperature reconstruction technique for the South Pole ice core, Kahle et al. (2021) find an interglacial site temperature warming of 7.15±0.68C between 19.5–22.5 and 0.5–2.5 ka. Our reconstruction yields a site temperature warming of 8.9±0.4C for the same interval.

Buizert et al. (2021) used a similar approach as Kahle et al. (2021) (but without the diffusion constraint) to reconstruct Antarctic temperatures for many of the sites we investigate. Their results for the temperature change since the Last Glacial Maximum (LGM) are in close agreement with ours for West Antarctica but are substantially smaller than ours in East Antarctica (and less than Kahle et al.2021, for the South Pole). We emphasize that the firn modeling approach cannot simultaneously satisfy all observational constraints, as discussed in Kahle et al. (2021), suggesting that the differences may lie in a still incomplete understanding of firn processes (see, e.g., Gkinis et al.2021). Alternatively, Buizert et al. (2021) suggest that changes in the inversion strength may explain the differences.

Figure 13Antarctic-wide stacks of reconstructed temperature histories. (a) Moisture source evaporation temperatures, T0, of all ice-core sites are shown in light blue, while the average of all records is shown in dark blue. (b) All reconstructions of ice-core site surface temperatures, Ts, are shown in light red, while the continent-wide average is shown in dark red. (c) Anomalies of site-averaged evaporation temperature (blue), condensation temperature (gold), and surface temperature (red) are shown with respect to the mean value of the most recent 2000 years. All records are interpolated to even 50-year spacing. Where possible records are on the synchronized Buizert et al. (2018) timescale.


We create a stack of each reconstructed temperature variable (the evaporation temperature T0, condensation temperature Tc, and surface temperature Ts) for all eight ice-core records (Fig. 13). We weight the records equally; we do not adjust for the spatial distribution of the cores or weight by area, latitude, or elevation.

In Fig. 13 one can see that the Antarctic-wide average surface temperature change during the last deglaciation was considerably larger than the concurrent temperature change in the mean moisture evaporation source. In fact, average deglacial change in Antarctic surface temperature was about 3 times as large as the changes in evaporation temperatures, while changes in condensation temperature were about twice as large as the evaporation temperature changes.

In Fig. 14 we show the pattern of glacial–interglacial temperature change across the Antarctic continent. The magnitude of warming since the Last Glacial Maximum is calculated as the temperature difference between the late Holocene (LH, 0–4 ka) and Last Glacial Maximum (19–23 ka) for comparison with other proxy reconstructions. There may be some uncertainty in the relative magnitudes of these changes owing to offsets in the individual timescales of each record. While the relative magnitudes of interglacial change depend on the exact time periods used in the differencing, the pattern of changes in surface temperature across the continent is robust. Disentangling the full suite of causes for these temperature changes is a topic of future work but may certainly include changes in the local energy balance, heat and moisture transport, and ice sheet topography.

Figure 14Spatial pattern of temperature change since the Last Glacial Maximum. (a) Warming between 19–23 and 0–4 ka for ice-core site surface temperatures (colored circles with black outline corresponding to different ice-core sites as shown in map inset, uncertainty in the temperature change shown as error bars), moisture source evaporation air temperature (colored circles, red outline and uncertainty), and moisture source sea surface temperature (red circles). Moisture source warming is plotted at the mean latitude of the modern moisture source distributions for each ice-core site, while the latitudinal extent of each moisture source is indicated by the relative histograms along the x axis. Sea surface temperature warming from the MARGO compilation of SST estimates from ocean sediment cores is shown as open black circles. (b) Ice-core site surface temperature changes and moisture source sea surface temperature changes are shown as large colored circles with black outline. MARGO compilation SST changes are shown as small colored circles. (c, d) Spatial pattern of temperature change from the multi-model mean PMIP3 simulations of the LGM and pre-industrial. (c) The multi-model mean for all grid points is shown as grey dots with the zonal mean in black. Estimates of ΔTs and ΔT0 from the ice-core reconstructions are shown as colored circles as in (a) for reference.

We find smaller glacial–interglacial temperature change for East Antarctic sites compared to previous reconstructions. Our results show that the surface temperatures of the lower, warmer areas of West Antarctica warmed significantly more than the higher, colder East Antarctic since the LGM. For example, the average surface temperature warming between the LGM and LH for the two lowest sites in our reconstruction, WDC and Siple Dome, is 11.6 C. The average warming at the two highest sites of Dome Fuji and Vostok, however, is significantly less and just 6.9 C or 59 % of that at the lower sites.

In Fig. 14 we plot the magnitude of warming since the LGM of Antarctic moisture source evaporation air temperatures for all ice-core records as a function of the mean latitude of the moisture source distribution for each site (based on a water-tagged general circulation model (GCM) simulations; see the Appendix). Additionally, we calculate the change in sea surface temperature (SST) during evaporation (red circles Fig. 14a) using the T0–SST relationship from our model for comparison to other SST proxy reconstructions. While plotted as points, note that these changes in moisture source temperature reflect the integrated warming over the moisture source area. Modern moisture source distributions for each site are indicated by relative histograms along the latitude axis in Fig. 14a (see Fig. 1 for further information), which are colored to correspond to the ice-core site. The moisture source points are plotted at the longitude of the respective ice-core sites in Fig. 14b, though in reality the moisture source distributions (MSDs) have a significant meridional extent (see extended data in Fig. 8, Buizert et al.2018) that is often asymmetrically to the west owing to the westerly winds. Changes in T0 for Antarctic moisture sources may reflect both warming SSTs at fixed locations and potential changes in the mean latitude of the moisture source distributions, for example due to changes atmospheric circulation, e.g., a meridional shift in the mean westerly winds (Markle et al.2017). Disentangling these two influences requires additional constraints and is beyond the scope of this study. While the ice-core T0 reconstructions have low spatial resolution owing to broad moisture source distributions, they benefit from the temporal resolution and precision of the ice-core age scales compared to other proxy records of temperature from the Southern Hemisphere midlatitudes.

It is clear from Fig. 14 that the Antarctic continent warmed 2 to 3 times as much as the Southern Hemisphere midlatitude moisture source areas since the Last Glacial Maximum. This result is in line with other paleoclimate reconstructions, as well as modeling of the pattern of polar amplification since the LGM (Masson-Delmotte et al.2006; Otto-Bliesner et al.2006). In particular, our estimates of moisture source region changes agree with completely independent estimates from the MARGO compilation of SST changes (MARGO Project Members2009, open circles Fig. 14). There appears to be some zonal asymmetry in the warming of Southern Ocean surface temperatures in both our moisture source reconstructions and the MARGO compilation. The waters around New Zealand and Australia that comprise the moisture source of Talos Dome appear to show the most warming since the LGM.

These patterns of Southern Hemisphere warming are also in reasonable agreement with modeling expectations, e.g., from the Paleoclimate Model Intercomparison Project (PMIP3; Braconnot et al.2012). The multi-model mean pattern of Southern Hemisphere polar amplification from PIMP3 simulations is shown in Fig. 14. There is broad similarity to our reconstructions, though there are important differences as well. The spread in temperature change about the zonal mean over both the Antarctic and ocean surface is similar between the model and the reconstructions. Our reconstructions show more warming in the ice-core moisture source areas equatorward of the polar front than the PMIP3 mean and less warming over the surface of West Antarctica. We note that the magnitude and pattern of modeled Antarctic surface warming is predominately a function of imposed changes in ice sheet surface elevation to PMIP3 experiments. The extreme warming seen in parts of West Antarctica in the PMIP3 model mean (e.g., >20C), which is outside the range found in our reconstructions, likely reflects unrealistically large ice sheet thickness changes prescribed in PIMP3. This assessment is consistent with the findings of Werner et al. (2018), who made an extensive analysis of GCM-modeled changes in Antarctic water isotopes and ice-core records. Understanding the full set of processes responsible for the reconstructed pattern of Southern Hemisphere polar amplification is a topic for future work.

6 Conclusions

Ice-core records of the stable isotopes of water provide detailed histories of Earth's climate. Both qualitative and quantitative interpretation of these records requires understanding the relationships between fractionation processes and environmental conditions.

Qualitatively, δ18O and δD are reliable indicators of the relative change in condensation temperature over a sufficiently long timescale. The assumption of a roughly linear relationship is generally justified, as shown in this study and previously. However, the linear definition of deuterium excess, dxs, is an unreliable indicator of relative evaporation site temperature change, particularly at East Antarctic sites with very depleted δ18O and δD values. In these cases, the logarithmic definition of the parameter, dln, is a more faithful qualitative proxy for evaporation temperature.

We can use models to make quantitative interpretations of water-isotope variability and to disentangle the combined influences of the source and site temperatures. To date, most water-isotope temperature inversions have assumed linear relationships (Kavanaugh and Cuffey2003; Vimeux et al.2002; Stenni et al.2011; Uemura et al.2012). However, as shown here, this assumption is flawed. Even in the simplified water-isotope models that underlie most temperature reconstructions, there are inherent nonlinearities in the isotope–temperature relationships. Ignoring these nonlinearities distorts reconstructed temperature variability. In the case of evaporation source temperature changes, these distortions may be a significant fraction of the total reconstructed variability.

There is a long-standing debate regarding the interpretation of “spatial” and “temporal” slopes in the water-isotope–temperature relationship (e.g., Jouzel et al.1997). These discussions are conceptually useful. However, while space and time are obvious coordinates through which to understand climate, they are not the most relevant for water-isotope fractionation. Neither space nor time can independently cause water to change phase and fractionate.

The fundamental dimension through which to understand water-isotope fractionation is temperature. In this study we use a relatively simple model to investigate the relationships of water isotopes in precipitation to temperature. While the distinction between temporal and spatial slopes is not directly addressed in this context, we are able to resolve the core question: is the water-isotope–temperature relationship fixed? It is not. But these slopes are fundamentally functions of temperature, not space or time.

Our nonlinear reconstruction technique allows for the estimation of absolute temperatures in the past, in addition to their variability, and is corroborated by independent temperature constraints. By taking into account the inherent nonlinearities of water-isotope fractionation we are better able to constrain evaporation source region changes. Our reconstructions reveal a spatial pattern of temperature change across the Antarctic continent in which West Antarctica warmed significantly more than East Antarctica since the Last Glacial Maximum. Further, our reconstructions provide insight into the spatial pattern of polar amplification, suggesting that the warming since the LGM in Antarctica was 2 to 3 times that in the midlatitudes.

Appendix A: Simple Water Isotope Model

The Simple Water Isotope Model (SWIM) is based on existing numerical Rayleigh-type distillation models (Merlivat and Jouzel1979; Jouzel and Merlivat1984; Ciais and Jouzel1994; Criss1999; Kavanaugh and Cuffey2003), though we make several important improvements and updates. We first describe the physical environmental aspects of the model and then the details of the fractionation scheme.

A1 Environmental trajectory

Our model considers moisture transported from evaporative sources down an atmospheric temperature gradient (i.e., from the midlatitudes toward the pole), driving condensation and fractionation. Our model operates in the dimension of temperature; we consider pseudo-adiabatic temperature pathways from an initial surface air temperature, T0, to a final condensation temperature, Tc, and discrete steps dT, as well as Euler numerics.

A1.1 Source region conditions

The moisture source surface air temperature (Ta), sea surface temperature (SST), and relative humidity (RH) influence the fractionation of vapor evaporating from the ocean. We use modern climatological correlations to find initial values of SST0 and RH0 given a specified initial air temperature, T0, using the 1980–2010 annual mean climatological fields from the NCEP/NCAR reanalysis project (Kalnay et al.1996) and the ERA-Interim reanalysis (Dee et al.2011). Correlations with surface air temperature are better defined than spatial correlations and give greater flexibility to the model for use in different climate states. Surface air temperature and SST are extremely well correlated in the reanalysis (Fig. A1) with a well-defined nearly linear relationship over most of the temperature range, except where the SSTs asymptote to the freezing point of seawater. The relationships are nearly identical between the NCEP and ERA reanalysis products.

Relative humidity gradients in the modern climate are fairly weak, though surface RH over the ocean is consistently higher at lower surface temperatures on climatological timescales. While variable on short timescales, RH appears to be largely invariant on timescales longer than interannual (Dai2006; Vimeux et al.2002). We find that the over-ocean surface relative humidity is systematically about 5 % higher in the NCEP reanalysis compared to the ERA reanalysis, though the relationship to surface temperature is similar.

Given a specified initial air temperature, T0, our model uses values of SST0 and RH0 based on fits to the modern climatology. We use three methods to calculate the climatological relationships over the interval −10Ta 28 C: a cubic spline with specified noise tolerance, the mean and median of SST and RH distributions within binned values of Ta, and high-order polynomial fits. All methods show effectively indistinguishable relationships in both reanalysis products. We calculate the uncertainty in these fits and test the model's sensitivity. Our base model uses the cubic spline method, which is least susceptible to edge effects and phase shifting, and it maintains a smooth first derivative.

We find some differences in the Ta-to-SST fit between the Northern and Southern Hemisphere for air temperatures between 5 and −15C, as seen in Fig. A1. In this study we will use the Southern Hemisphere fit for Ta to SST. We find no major hemispheric differences for the Ta-to-RH fit and find little impact of zonal asymmetry on either fit. We find relatively small differences in the fit between Ta and SST0 for different seasons and somewhat larger seasonal changes in the Ta and RH0 relationship. We use the annual average fits hereafter and test the sensitivity of the water-isotope values of evaporation to these seasonal differences.

Figure A1Climatological correlations between Ta, SST, and RH. (a) Annual mean climatological surface air temperature, Ta, and sea surface temperature, SST, from the NCEP/NCAR reanalysis (Kalnay et al.1996). Light blue dots are Northern Hemisphere (NH) grid points, while dark blue dots are Southern Hemisphere (SH) grid points. The polynomial fit for the SH is in red. The error estimate for the fit, σerr in orange, is the standard deviation of the misfit in the model. The 1:1 line is shown in black. (b) Same as on the left but for the surface air temperature over ocean, Ta, and surface relative humidity (RH) over oceans.


The normalized relative humidity, RHn, is critical to kinetic fractionation of water isotopes during evaporation (Merlivat and Jouzel1979; Risi et al.2010) and depends on the three variables above:

(A1) RH n = RH × e s ( T a ) e s ( SST ) ,

where es (Ta) and es (SST) are the saturated vapor pressures of air at the surface air temperature and at the sea surface temperature, respectively.

A1.2 Transport

After evaporation at initial air temperature, T0, and specified surface pressure, P0, moisture is transported toward the pole in isolation, cooling and condensing along the way. The air parcel is cooled pseudo-adiabatically, defining a pressure trajectory, P as a function of temperature, T (Fig. 2). As the air parcel cools, moisture above saturation is removed and the latent heat released during the phase change keeps the air parcel warmer than in an otherwise equivalent isobaric pathway. Following the pseudo-adiabatic assumption, we consider no other heat sources to the air parcel, and moisture is removed immediately after condensation. Below we investigate in the influence of air parcel mixing on our model framework, which is a relaxation of the adiabatic assumption.

We calculate a pseudo-adiabat following the iterative routine described in Bakhshaii and Stull (2013) but taking into account the saturated vapor pressures of both ice and liquid water condensate. The temperature-dependent saturated vapor pressures of ice and liquid water (Murray1966; Murphy and Koop2005), together with air pressure P(T), define saturated mixing ratios for ice and liquid water,

(A2) r s = R d R wv × e s P - e s ,

as functions of T, where RdRwv is the ratio of gas constants of dry air and water vapor.

We consider air parcels with mixed ice and liquid condensate (Ciais and Jouzel1994), in which the ice fraction smoothly increases as temperatures decreases below freezing. Many models, including isotope-enabled GCMs, approximate the temperature dependence of cloud ice–liquid fraction as piecewise linear functions (Hu et al.2010), while others use smoothly varying error integrals (Ciais and Jouzel1994; Kavanaugh and Cuffey2003). We use temperature-dependent functions for the cloud ice fraction derived from satellite observations (Hu et al.2010) over the Southern Ocean and over land snow and ice surfaces, which preserve significantly more liquid water at colder temperatures than previous parameterizations (e.g., Kavanaugh and Cuffey2003).

The specific heat at constant pressure, cp, latent heat, L, saturated vapor pressure, es, and saturated mixing ratio, rs, are temperature-dependent and calculated for the liquid and ice phases individually. Effective values for the parcel as a whole are calculated from the mixing fractions of each phase (Kavanaugh and Cuffey2003). For example, rs(eff)=rs(ice)F(ice)+rs(liq)F(liq), where rs(ice) and rs(liq) are the saturated mixing ratios of ice and liquid, respectively, and F(ice) and F(liq) are the (temperature-dependent) fractions of each phase of condensate.

Moisture is removed along the temperature pathway owing to the temperature-dependent changes in the saturated mixing ratio, -dqdT=drs(eff)dT. This is a simplified view of large-scale precipitation commonly used in similar models (e.g., Markle et al.2018). We do not consider reevaporation of falling precipitation. There are several reasonable choices in the implementation of our simplified view of moisture removal once air parcels are cooled from initial relative humidity to saturation. The instantaneous moisture removal process may leave the air parcel at saturation, at some specified level below saturation (e.g., RH = 90 %), or at the air parcel's initial level below saturation, in which case relative humidity is constant along the path. We test the sensitivity of our model to these assumptions below using constant relative humidity as our default.

At very cold temperatures moisture is removed not at saturation but at a specified level of supersaturation. The presence of both ice and liquid condensate in the cloud dictates a supersaturation of vapor over ice due to the difference in liquid and ice vapor pressures (Jouzel and Merlivat1984). A paucity of condensation nuclei may lead to further supersaturation at very cold temperatures (Tegen and Fung1994). The total supersaturation is parameterized here to depend on temperature (discussed in detail in Sect. A2.3).

A2 Isotope fractionation

In this section we outline the water-isotope fractionation scheme used in SWIM. We model equilibrium and the kinetic fractionation of the 2H1H, 18O16O, and 17O16O ratios in water. We use conventional notation in which R is the number ratio of heavy to light isotopes of a species, for example DR=2H1H and 18R=18O16O.

The fractionation factor is the ratio of R values between phases. For example, the fractionation factor between liquid and vapor phases for δ18O is

(A3) 18 α l - v = 18 O 16 O liquid 18 O 16 O vapor = 18 R l 18 R v .

We use the empirically determined, temperature-dependent equilibrium fractionation factors between liquid and vapor, 18αeq(l−v) and Dαeq(l−v), as well as those between vapor and ice, 18αeq(i−v) and Dαeq(i−v) (Majoube1970, 1971; Merlivat and Nief1967; Criss1999), with updates for the ice–vapor equilibrium fractionation factor found by Lamb et al. (2017).

A2.1 Evaporation from the ocean

The isotopic values of vapor evaporating from the ocean are determined, in part, by the isotopic values of the seawater. By definition, globally averaged seawater has δ values near 0 ‰. However, the δ18O of seawater (δ18Osw) in the regions of Antarctic moisture sources is more depleted than average, with a mean around −0.3 ‰, (Schmidt et al.1999). We use the observed correlation between δ18Osw and δDsw from a compilation of global seawater measurements (Schmidt et al.1999) to find an initial δDsw given the specified initial δ18Osw (which changes with mean climate) and investigate the sensitivity of the model to these initial conditions. We assume a 17Oxs of seawater equal to 0, where 17Oxs=δ17O-0.528×δ18O, which is typically reported in per meg.

The atmosphere above the global oceans is not at saturation on average, with relative humidity typically around 80 % (Hartmann2015). Because of this steady-state disequilibrium, significant kinetic fractionation occurs during evaporation from the ocean. Kinetic fractionation depends both on the relative humidity and the wind speed at the air–ocean interface during evaporation (Merlivat and Jouzel1979). The effective fractionation factor associated with diffusion and turbulence is

(A4) α diff = D D * n ,

where D and D* are the diffusivities of the light and heavy isotopes, respectively (Merlivat and Jouzel1979). The exponent n ranges from 0 to 1 and relates to the wind regime and speed and the ratio of turbulent to molecular diffusion. For the diffusive fractionation between H218O and H216O during initial evaporation, the fractionation factor 18αdiff equals 1.0 for pure turbulence and 1.0028 for pure molecular diffusion (Merlivat and Jouzel1979; Barkan and Luz2007).

Following Kavanaugh and Cuffey (2003), we do not explicitly consider surface wind speeds. Instead we use an empirical approach and the results of Uemura et al. (2008, 2010) for 18αdiff, who estimate the parameter based on measurements of δD, δ18O, and δ17O in vapor above the Southern Ocean. Uemura et al. (2010) find a value of 18αdiff=1.007±0.0013 and 1.008±0.0018 when optimizing for observations of dxs and 17Oxs of vapor, respectively. These results are within uncertainty of each other and of independent analysis by Pfahl and Wernli (2009), which found a value of 1.0076. Using a compilation of vapor measurements (including Uemura et al.2008, 2010; Liu et al.2014; Kurita et al.2016; Benetti et al.2017), we find that 18αdiff=1.009 leads to a good match between modeled and observed values of both dxs and 17Oxs of vapor when using the observed values of T0, SST, and RH at the time of the vapor measurements. We investigate the sensitivity of the model to this parameter below.

The diffusive fractionation factor between hydrogen and deuterium, Dαdiff, may be determined experimentally by measuring the ratio of diffusive fractionation factors (Merlivat1978; Luz et al.2009).

(A5) ϕ diff = D α diff - 1 18 α diff - 1

Merlivat (1978) found a mean value for ϕdiff of 0.88 based on laboratory evaporation studies, and Luz et al. (2009) found that the value of ϕdiff depends on the evaporation temperature, ranging between 0.73 and 1.06 for temperatures between 10 and 69.5 C. We use a piecewise linear function based on the results of Luz et al. (2009) to relate ϕdiff and evaporation temperature and thus Dαdiff to 18αdiff. For evaporation temperatures colder than the experimental range of Luz et al. (2009) (<10C), we use the measured value at 10 C (ϕdiff=1.06). The differences in model results for evaporation with a temperature-dependent ϕdiff and a constant ϕdiff=0.88 are small (<1 ‰ for initial δD of vapor).

For the fractionation of 17O16O we use the following relationships, which are backed by both theory and empirical observation (Barkan and Luz2005, 2007): 17αeq=18αeq0.529 for vapor and liquid in equilibrium and 17αdiff=18αdiff0.518 for vapor diffusion.

An alternative to the empirical approach for diffusive fractionation outlined above is to use new results from kinetic molecular theory (Hellmann and Harvey2020), which lead to physically constrained temperature dependence of the diffusivity ratios for the isotopologues. These results lead to a temperature dependence of ϕdiff that is different than the experimental results of Luz et al. (2009), being less variable and closer to the 0.88 value of Merlivat (1978). Using this formulation in the evaporation scheme requires us to choose a representative value of n (see Eq. A4) for the wind turbulence regime to match observations. A value of n between 0.22 and 0.32 leads to a range of results equivalent to that resulting from 18αdiff between 1.006 and 1.010 and the empirically based scheme described above. As before, this temperature dependence leads to very minor differences for the δ values of vapor as a function of temperature.

The relationship between the initial R value of the vapor and the ocean due to kinetic fractionation depends on the normalized relative humidity during evaporation, RHn, and the equilibrium and diffusive fractionation factors, αeq(l−v) and αdiff. Following Criss (1999) and Luz et al. (2009),

(A6) α evap = R o R e = α eq α diff ( 1 - RH n ) 1 - α eq RH n R v R o ,

where Ro and Rv are the isotopic ratios of the ocean water and the water vapor in the atmospheric boundary layer, respectively. Re is the ratio of the evaporate, the net vapor lost to the atmosphere, which is a quantity that is not directly measurable (Criss1999).

If we assume that the only source of vapor to the boundary layer is the local evaporate, we may equate Rv and Re and solve Eq. (A6) for Rv (Merlivat and Jouzel1979; Criss1999; Risi et al.2010):

(A7) R v = R o α eq × α diff + RH n 1 - α diff .

The modeled isotopic composition of vapor evaporated from the ocean is shown in Fig. A2. This “local” closure assumption is within the range of observations of water isotopes in vapor over the Southern Ocean and elsewhere (Uemura et al.2008, 2010; Liu et al.2014; Kurita et al.2016; Benetti et al.2017).

However, the validity of the local closure assumption under certain conditions is in question (Uemura et al.2010; Risi et al.2010). In addition to moisture from the ocean surface, the boundary layer may receive moisture from advection, convection, subsidence, and reevaporation of precipitation. Risi et al. (2010) explored this issue extensively using a model that takes into account these other sources of moisture. They show that the local closure assumption leads to vapor that is too enriched in both δD and δ18O and too low in dxs, and these offsets are a function of environmental conditions (Risi et al.2010).

Figure A2Modeled isotopic values of evaporation compared to Southern Ocean vapor measurements. (a) Modeled δ18O versus SST using the “local” closure assumption (red lines) and “global” closure assumption (blue lines). We use two different reanalysis datasets for the SST and RH climatology: NCEP/NCAR (solid lines) and ERA-Interim (dashed lines). Black dots are discrete vapor measurements from the Southern Ocean made by Uemura et al. (2008) (U08), while grey, blue, red, and yellow dots are continuous ship-based Southern Ocean vapor measurements made by Liu et al. (2014) (L14), Kurita et al. (2016) (JARE55), and Benetti et al. (2017) (STRASSE, RARA AVIS), respectively. (b) Same as (a) but for modeled δD versus SST. (c) Same as (a) but for modeled dxs versus SST. The vertical tails at low SST in panels (a) and (b) result from SSTs asymptoting to the freezing point of seawater while air temperatures may continue to decrease.


We investigate the influence of the closure assumption in a few ways. First, we examine closure globally rather than locally. The mean ocean has δ values of about 0 ‰, and global average precipitation has δ18O=−4.5 ‰ and δD=-26.7 ‰ (Craig and Gordon1965). Considering the global average steady state in which the ocean is the ultimate vapor source for precipitation (Merlivat and Jouzel1979; Criss1999) , the delta values of precipitation must reflect the net loss by evaporate from the ocean. Thus, globally, the ratio RoRe is 1.0267 for D and 1.0045 for 18O. Instead of equating Re and Rv locally, we define a global αevap=RoReglobal and solve for Rv. Substitution into Eq. (A6) and rearranging leads to

(A8) R v = R o 1 - α eq α diff 1 - RH n α evap α eq RH n - 1 .

Though an obviously blunt approach for determining local evaporation, this global closure assumption is the limit for a globally mixed atmosphere. Modeled evaporation using both closure assumptions is compared to isotopic measurements of Southern Ocean vapor in Fig. A2.

Next, we consider specific mixing of evaporative conditions instead of the generalized globally mixed case above. Rather than the local ocean being the only source of vapor, we can consider a simple scenario in which the isotopic composition of the boundary layer, Rv, is comprised of both local evaporate, Re, and vapor evaporated at some distal location and advected to the site, Rv=(1-θ)Re+θRdistal, where θ is the fraction of nonlocal vapor with composition Rdistal. The local evaporative conditions are defined by T0, while the distal conditions may be either warmer or colder than T0. Vapor is evaporated under those distal conditions (using a local closure assumption) and advected without fractionation before mixing with the local vapor of composition Re. Figure A3 shows the isotopic composition of vapor over a range of T0, with a full range of contributions from both a 5 C warmer moisture source (T0+5C, red) and a 5 C colder moisture source (T0−5C, blue). This range of mixing leads to a spread of delta values of the initial vapor around the simple local closure assumption, though the difference is generally less than that between the local and global closure assumptions.

Figure A3Relationship between isotopic composition of vapor and SST under different mixing scenarios. Positive values on the color scale reflect the fraction of moisture from a 5 C warmer moisture source mixed with the local moisture source (50=50 % moisture from local source + 50 % moisture from 5 C warmer-than-local moisture source), while negative values reflect the fraction of moisture from a 5 C colder moisture source mixed with the local moisture source (-50=50 % moisture from local source +50 % moisture from 5 C colder-than-local moisture source). Model results use NCEP/NCAR reanalysis for SST and RH climatology. Unmixed model results for a local closure assumption are shown as a black dashed line and a global closure assumption as a solid black line.


The isotopic values of vapor produced by any of these closure assumptions (local, mixed, or global) are within the range of mean values of the observational data and show similar relationships to local environmental conditions like SST and RH. These closure assumptions represent the bounds of a well-mixed and unmixed atmosphere or something in between. The amount of mixing in the boundary layer could change with location and with climate mean state. Rather than tying our model to one closure assumption, we view mixing at evaporation as an inherent uncertainty.

Figure A4Comparison of modeled and observed isotope excess parameters and relationship to source region conditions. (a) Observed dxs and SST relationship in Southern Ocean vapor from Uemura et al. (2008) (black dots, U08) and Liu et al. (2014) (grey dots, L14). SWIM results for evaporation under SST and RH conditions observed coincident with vapor measurements of Uemura et al. (2008) (cyan dots, U08 model) and Liu et al. (2014) (purple dots, L14 model). (b) Same as (a) but for modeled and observed dxs-to-RH relationship from observations of Uemura et al. (2008) The Liu et al. (2014) observations and model show a similar trend and are omitted for visual clarity. (c) Observed 17Oxs and SST relationship in Southern Ocean vapor from Uemura et al. (2010) (black dots, U10) and SWIM results run under observed sea surface conditions (cyan dots, U10 model). (d) Same as (c) but for observed and modeled 17Oxs and RH relationship in Southern Ocean vapor.


It is important to note that in Fig. A2 we use climatological correlations between T0, SST, and RH, while the observational data represent far shorter time intervals, mostly from one season. When using the observed values of T0, SST, and RH at the time of the observational measurements (Uemura et al.2008; Liu et al.2014), we are able to capture the complex variability in the isotopic values of the vapor on those given days. For example, in Fig. A4 we compare modeled dxs of vapor and modeled 17Oxs of vapor to Southern Ocean vapor observations using the observed environmental conditions at the time of the vapor measurements. The modeled relationships between dxs and 17Oxs with SST0 and RH0 are in good agreement with observations.

We examine the sensitivity of initial evaporation to several model parameters discussed above in Fig. A5. In all cases the modeled sensitivity to these parameterizations and uncertainties is relatively small compared to the natural variability in observations of isotopic vapor compositions. The choice of reanalysis product used to derive the climatological relationships between T0, SST0, and RH0, as well as the uncertainty in those fits, has relatively small effects on the results of evaporation (Figs. A2, A5a–f). We also show the influence of the initial δ18Osw of the ocean (Fig. A5g–i) as well as the value of αdiff (Fig. A5j–l).

Figure A5(a–f) Sensitivity of modeled δ18O, δD, and dxs of vapor to uncertainty in the reanalysis-based fits between climatological Ta, SST, and RH in the NCEP/NCAR reanalysis. Panels (a) and (b) show the modeled isotope vapor relationship to uncertainty in the climatological Ta-to-SST relationship. Red lines show the model run using the central estimate of the fit, and orange lines show the spread expected with ±σerr of the fit as shown in Fig. A1. Panels (c) and (d) are the same as (a) and (b) but showing the central estimate (dark blue) and spread associated with ±σerr (light blue) in the climatological Ta-to-RH relationship shown in Fig. A1. (g–i) Sensitivity of modeled isotope values of vapor to the δ18Osw of seawater. Values of δ18Osw from −0.5 ‰ to 0.5 ‰ are specified, representing most of the global variance in δ18Osw. Values of δDsw are determined based on correlations of δ18Osw and δDsw from observations (Schmidt et al.1999). (j–l) Sensitivity of modeled isotope values of vapor to a range of 18αdiff values from 1.007 to 1.010. Also shown are results using the first-principles-based diffusivities of Hellmann and Harvey (2020) (HH20, black, using n=0.27 in Eq. A4). Local closure is assumed in all panels.


Figure A6Relationship between observed and model dxs (a) and 17Oxs (b) of vapor over the Southern Ocean. Observed vapor values are from Uemura et al. (2008, 2010) (large dots) and Liu et al. (2014) (small dots). Modeled values use the reported T0, SST0, and RH0 from the observations and four different values of 18αdiff (shown in color of dots). The 1:1 line is shown in black.


The direct comparisons of observed and modeled vapor composition using observed T0, SST0, and RH0 at the time of the vapor measurements (Fig. A4) suggest that a single effective αdiff may not fully capture the kinetic effects across the range of surface conditions. While αdiff=1.009 leads to a good fit between observed and modeled deuterium excess for much of the range of surface conditions, there is a small persistent misfit for surface temperatures between 20 and 27 C, where a smaller αdiff is suggested (Fig. A6). While it is possible to implement a temperature-dependent αdiff to reduce this misfit, we prefer a fixed αdiff to avoid over fitting a relatively small dataset in the absence of further evidence or physical reasoning. While we do not consider temperature dependence of αdiff, we do consider a range of αdiff as an inherent uncertainty in our model and account for this in the uncertainty analysis of our temperature reconstructions as discussed in Sect. A9.

We note that it is of course also possible that other factors, rather than temperature dependence of αdiff, could account for the apparent misfit, such as a difference between the ship-measured RH and Ta and that felt at the water's surface or specific mixing of nonlocal moisture in the boundary layer on the days of the ship-based measurements. Given the magnitude of the misfit, it is also possible that spatial or temporal variability in δ18Osw and δDsw could account for the misfit (Schmidt et al.1999).

A2.2 Distillation

We next discuss the distillation of water isotopes during transport. As an air parcel cools, water condenses, fractionates, and is removed as precipitation. The essential differential equation for Rayleigh distillation (Rayleigh1902; Dansgaard1964; Criss1999) is

(A9) d ln ( R ) d ln ( f ) = α - 1

where f is the fraction of initial water vapor remaining in the air parcel. The amount of moisture at any temperature along the pathway is found by integrating the changes in the saturated mixing ratio rs owing to pseudo-adiabatic cooling from the source (Dansgaard1964) . Thus,

(A10) f = q q 0 = r s r s 0 .

In general, condensation occurs in the model at saturation, and thus the temperature-dependent equilibrium fractionation factor αeq is used in Eq. (A9). However, at cold conditions there may be supersaturation of vapor over ice, leading to additional kinetic fractionation. Following previous models (Jouzel and Merlivat1984), the total fractionation αtot factor is αtot=αeqαk. Equation (A9) thus becomes

(A11) d ln ( R ) = ( α tot - 1 ) d ln ( f ) .

The kinetic fractionation factor, αk, depends on the supersaturation of vapor over ice, Si:

(A12) α k = S i α eq × D D * S i - 1 + 1 .

Following Jouzel and Merlivat (1984), we use the ratio of diffusivities for oxygen isotopes D16D18=1.0285 during moisture transport, representative of pure molecular diffusion and ignoring the negligible ventilation effect. Likewise we use D1D2=1.0251 for the ratio of diffusivities of hydrogen isotopes. These values imply a constant ϕdiff during transport equal to 0.88 (Jouzel and Merlivat1984) rather than the temperature-dependent ϕdiff used in the evaporation scheme. We prefer this value for simplicity, consistency with earlier work, and lack of experimental measurements of ϕdiff at the colder temperatures experienced during transport (Luz et al.2009). However, to investigate the effect of assuming a constant ϕdiff, we also use a configuration of the model with the temperature-dependent diffusivity ratios of Hellmann and Harvey (2020) based on molecular kinetic theory. We can achieve essentially identical results using this physically based temperature dependence of the diffusivity ratios and assuming no temperature dependence, with exceedingly minor changes to the supersaturation function (see Sects. A2.3 and A4); a change of b in the supersaturation function of less than 0.00015 is sufficient to negate any difference between these assumptions. We thus consider the effect of temperature dependence in these diffusivities to be subsumed within the uncertainty associated with tuning the supersaturation function (see Sects. A4 and A9)

In the mixed-phase portion of the transport pathway, the effective fractionation factors are determined by the mixing fractions of ice and liquid condensate. Following Kavanaugh and Cuffey (2003),

(A13) α eff = α tot ( l - v ) F ( liq ) + α tot ( i - v ) F ( ice ) .

The temperature dependence of the ice and liquid fraction is shown in Fig. A7a and based on satellite observations (Hu et al.2010).

A2.3 Supersaturation

The supersaturation of vapor over ice is a critical parameterization in water-isotope distillation models. The true relationship of supersaturation to environmental conditions is the result of complex cloud microphysics (Hong et al.2004). Because of its strong influence on water-isotope fractionation and the uncertainty in the underlying physics, the supersaturation is often parameterized to depend on temperature and tuned to fit water-isotope models to observations (Jouzel and Merlivat1984; Kavanaugh and Cuffey2003; Schoenemann et al.2014). Jouzel and Merlivat (1984) parameterized the supersaturation as a function of temperature and note that available water-isotope data could not distinguish among possible functional forms of the parameterization (e.g., linear, exponential). Their linear parameterization has been used extensively in water-isotope models since:

(A14) S i = a - b × T ,

where a and b are tuned to fit observational data.

It is important to note here that the prescribed mixing of liquid and ice in the cloud implies a supersaturation of vapor over ice that follows the blue curve shown in Fig. A7b, which is inconsistent with the supersaturation driving kinetic fractionation as prescribed in Eq. (A14). The presence of both liquid and ice phases in a cloud is not the only potential source of supersaturation. The lack of condensation nuclei, for example, allows supersaturation to remain high in cold, ice-only conditions (Hong et al.2004), rather than returning to unity as the cloud becomes entirely ice-phase. It is common for water-isotope models, even those in some GCMs (e.g., Schoenemann et al.2014), to have multiple variables equivalent to supersaturation in different aspects of the same model, such as the isotope fractionation and precipitation schemes, which may not be self-consistent. Because the environmental supersaturation experienced by the air parcel is related to the relationship between temperature and moisture removal (that is dln(f)dT) and the supersaturation driving kinetic fractionation relates temperature to δ values (largely through Eq. A12), an inconsistency in the model's view of supersaturation can influence the modeled water-isotope–temperature relationship in unphysical ways.

To resolve this physical inconsistency, precipitation only occurs in SWIM when the parcel reaches the prescribed supersaturation by dictating an effective saturated mixing ratio of the air parcel, in which rs(eff)=rs(ice)Si. Ensuring consistent supersaturation across the model leads to a smoother relationship between temperature and the δ values of precipitation. This is in contrast to rather complex curvature in the temperature–water-isotope relationship that results if inconsistent relationships between saturation and temperature are used in the precipitation and fractionation schemes, which is generally incompatible with observations. In line with previous work (Jouzel and Merlivat1984), we find that using only the supersaturation implied by the mixing of ice and liquid, across all aspects of the model, results in a relationship between δ18O and δD irreconcilable with observations. Were the air parcel to return to non-supersaturated conditions in the ice-only portion of the cooling pathway, the simultaneous transition to equilibrium-only fractionation would drive a slope of δ18OδD that is not compatible with measured values in Antarctic precipitation. This gives additional credence to sustained supersaturation at cold temperatures.

Figure A7Supersaturation of vapor over ice. (a) Fractions of ice (cyan) and liquid (blue) condensate as a function of temperature. Curves are derived from satellite-based measurements (Hu et al.2010). (b) Supersaturation as a function of temperature. The blue curve shows supersaturation based solely on the saturated vapor pressures of ice and liquid, the mixing fractions based on the curves shown in panel (a), and the pseudo-adiabatic assumption (ILS). The red curve shows the linear parameterization of supersaturation (linear) used in the model, Si=1-0.00525C-1×T. Note that the atmosphere is subsaturated with respect to ice at temperatures above 0 C but that our pathways do not include ice at these temperatures.


A3 Application to Antarctic ice-core sites

A3.1 Moisture source distributions

The atmospheric circulation transports moisture poleward of ≈30 S (Fig. 1a). The mean evaporation latitude of moisture that precipitates at any given site can be estimated from moisture-tagged GCM experiments (Markle et al.2017). The difference between the mean latitudes of moisture evaporation and precipitation steadily increases between the subtropics and the pole (Fig. 1a and b). The mean evaporative latitude of moisture that precipitates in Antarctica ranges between 44 and 50 S (Fig. 1c). The surface elevation of the ice sheet is a strong predictor of the mean latitude of precipitation, with higher-elevation sites having more equatorward moisture sources (Fig. 1d) due to topographic isolation (Sodemann and Stohl2009; Bailey et al.2019). There are some notable asymmetries in this general pattern. The large embayments are areas of comparatively high-latitude mean moisture source, such as the Victoria Land coast in the Ross Sea region.

The mean moisture source latitude is, however, not the full story. The moisture reaching any Antarctic site does not originate from just a single mean source latitude or follow a single temperature pathway. The contribution of moisture evaporated from different latitudes to the final precipitation at a site defines a moisture source distribution (Markle et al.2017), which reflects the combination of the spatial pattern of evaporation, cumulative rainout, and the influence of atmospheric circulation. Here we diagnose annual mean moisture source distributions (MSDs) as a function of latitude from a moisture-tagged run of the Community Atmosphere Model (CAM) for East and West Antarctic sites (details are given in Markle et al.2017), as shown in Fig. 1e. Moisture source distributions derived from other methods like trajectory modeling are similar (e.g., Sodemann and Stohl2009; Markle et al.2012; Buizert et al.2018). These MSDs dictate the influence of evaporation source conditions (Ta, RH, SST) on moisture reaching any Antarctic site. There are some zonal asymmetries in surface conditions over the Southern Hemisphere oceans, but the strong latitudinal gradients are the largest source of spatial variance in these conditions at climatological timescales.

While the mean latitudes of moisture sources vary between Antarctic sites, largely as a function of site elevation, Antarctic MSDs are not fundamentally distinct in latitude, but rather span broadly overlapping swaths of the Southern Hemisphere from the Antarctic coast to the tropics (Fig. 1e). The difference in weighted mean moisture source latitude between Antarctic ice-core sites is less than 10 of latitude, while the moisture source distribution for any one site spans over 40 of latitude. Local evaporation is a small contribution to the moisture precipitating at Antarctic sites. On average moisture is transported more than 20 of latitude to reach Antarctica.

Given the broad range of evaporative conditions that contribute to moisture precipitating at an ice-core site, what is the meaning of the Tsource that can be reconstructed from water-isotope records? It is the moisture-weighted evaporative temperature, determined by the convolution of the spatial pattern of the MSD and the underlying surface temperatures (Fig. A28, Markle et al.2017). Both surface temperatures at fixed locations and the pattern of the MSD can change independently. The water-isotope records alone do not allow the disentanglement of these two patterns, which may have different temporal evolution (Markle et al.2017).

To understand the moisture transport and water-isotope distillation to Antarctic sites it is important to consider evaporation from the range of conditions comprising the moisture source distribution. We thus use an ensemble of temperature pathways for Antarctic precipitation defined by a range of Antarctic condensation temperatures as well as the broad range of evaporation temperatures underlying the Antarctic moisture source distributions. The means of these distributions vary across the continent.

A3.2 Condensation site conditions

During transport, moisture is cooled from initial surface air temperature at evaporation to subsequent condensation temperatures. The condensation temperature is not the same as the surface temperature where that precipitation falls. Indeed, there is a vertical and temporal distribution of condensation contributing to precipitation that falls at any point on the surface, analogous to the horizontal and temporal distribution of evaporation contributing the moisture ultimately transported to any precipitation site. What is the meaning of Tc reconstructed from ice-core records? It is the vertical profile of temperature weighted by the vertical profile of condensation that yields net accumulation to a site. The weighted condensation temperature has a distinct relationship to the surface temperature across the globe.

Antarctica has strong climatological inversions such that temperature aloft is often warmer than the surface (Connolley1996). Masson-Delmotte et al. (2008) review the relationship between the condensation temperature and the surface temperature (Ts) over Antarctica and compare the surface temperature to the weighted annual mean condensation temperature in both ERA-40 reanalysis (1980–2002) and MAR, a high-resolution mesoscale model forced by ERA-40 (see Fig. 8, Masson-Delmotte et al.2008). In both models the upper bound of the Antarctic condensation temperature appears to be set by the peak inversion temperature, though condensation temperatures are on average colder than the peak inversion temperature (meaning simply that condensation occurs at a range of temperature up to the peak inversion temperature). Masson-Delmotte et al. (2008) calculate a best fit of the surface to condensation temperature slope of 0.65 C C−1 in the ERA-40 data, consistent with previous work that found a slope of 0.67 C C−1 (Connolley1996; Jouzel and Merlivat1984). The spread of condensation temperatures in the higher-resolution MAR model suggests colder condensation temperatures than in the lower-resolution reanalysis (Masson-Delmotte et al.2008). The strength of the Antarctic inversion diminishes with increasing surface temperature (Connolley1996), and relatively warm Antarctic surface temperatures (e.g., >-20C) are associated with condensation temperatures colder than the surface temperature (Masson-Delmotte et al.2008).

We analyze the relationship between surface temperature and condensation temperature in monthly MERRA-2 reanalysis from 2008 through 2017 (Gelaro et al.2017). We show the climatological zonal-mean vertical profile of air temperature, the sum of the convective and large-scale precipitation source production rate, and the rate of reevaporation and sublimation of precipitation in Fig. A8. The relationship between the climatological weighted condensation temperature and the surface air temperature at every grid point is shown in Fig. A9. Note that this calculation accounts for the seasonality of precipitation throughout the atmospheric column, as well as the reevaporation and sublimation of falling precipitation. Ignoring reevaporation and sublimation leads to qualitatively similar results.

The primary take-away is that the MERRA-2 data show a generally linear relationship between condensation and surface temperature for typical Antarctic surface temperatures. That relationship, however, is not linear at warmer surface temperatures. Indeed, even at surface temperatures below zero, the relationship is not strictly linear, but rather steepens with decreasing temperature. The relationship between the surface air temperature and the weighted condensation temperature (for surface temperatures below −10C) has an average slope between 0.61 and 0.64 C C−1 depending on whether one accounts for reevaporation and whether the comparison is between the surface or 2 m air temperature. Note that this slope is weighted toward the surface temperature of regions comprising more model grid points. Further, the slope clearly steepens with decreasing temperature, reaching ≈0.710.75CC at the very coldest Antarctic surface temperatures. Given the uneven distribution of grid points in temperature space, it is difficult to estimate the robustness of this steepening of slope.

Figure A8MERRA-2 reanalysis data for 2008–2017. (a) Annual mean zonal-mean air temperature as a function of pressure and latitude. (b) Annual mean zonal-mean precipitation source (kilograms of water per kilogram of air per second). (c) Annual mean zonal-mean reevaporation and sublimation of falling precipitation in the same units and color scale as panel (b).


Figure A9Weighted condensation temperature as a function of surface temperature. (a) MERRA-2 reanalysis data for 2008–2017. Seasonally and vertically weighted condensation temperature (accounting for reevaporation and condensation) for every grid point in the Southern Hemisphere against 2 m air temperature (black) and surface temperature (dark grey). The same relationship with 2 m air temperature is shown for the Northern Hemisphere as light grey dots. (b) The reported surface temperature and the SWIM-modeled condensation temperature using pairs of δ18O and δD in the Masson-Delmotte et al. (2008) database (MD08, blue dots) from both Southern Hemisphere sites (black dots) and Northern Hemisphere sites (light grey dots) in the GNIP database, as well as the average of the top 50 years of sites from several deep ice-core sites (red dots). Solid black lines are 1:1.


Using our nonlinear temperature reconstruction method, we model the condensation temperature for every pair of δ18O and δD samples in the MD08 and GNIP datasets (Masson-Delmotte et al.2008; IAEA2001) that also have a reported mean surface temperature. We compare the relationship between the modeled condensation temperature and the reported surface temperature in Fig. A9. The pattern of this reconstructed relationship is remarkably consistent with that found in the MERRA-2 dataset, even to some extent at warm surface temperatures. This is despite the potential for offsets between the reported surface temperature and that at the time of the precipitation for the MD08 and GNIP precipitation samples, differing time intervals, potential moisture biases in the column in MERRA-2 (Bosilovich et al.2017), and the lack of processes in our isotope distillation model that should be important, for example, in tropical convection or that might alter Antarctic precipitation after deposition.

Examining all modeled condensation temperatures for samples in the MD08 and GNIP datasets with reported surface temperatures below 15 C, we find slopes between 0.62 and 0.67 C C−1. For just the Antarctic precipitation samples in the MD08 dataset we find a best-fit slope between the reported surface temperature and our modeled condensation temperature of 0.67–0.69 C C−1 (Fig. A9), depending on the model assumptions (0.69 C C−1 under our base assumptions) and whether below-freezing source evaporation is included (see below), in good agreement with previous Antarctic observational studies (Connolley1996; Jouzel and Merlivat1984). These slopes sit well within the range found in the MERRA-2 data. We also reconstruct the condensation temperatures for the topmost samples from several deep ice cores and compare those to the reported annual average temperatures for those sites (Fig. A9). We find a best-fit slope between 0.68 and 0.70 C C−1, depending on whether we average samples from the last 50 or 100 years, though only five points describe these lines.

Based on the above results we use the equation Tc=0.69CCTs-8.2C as our base estimate to reconstruct Antarctic surface temperatures,; however, we consider an uncertainty of ±0.02CC in the slope. Our base estimate leads to good agreement with the observed relationship between global δ18O and surface temperature (Fig. A10). We show our temperature-dependent isotope slopes in Fig. A11.

Figure A10(a) Observed relationship between δ18O of precipitation and reported surface temperature in the Masson-Delmotte et al. (2008) database (MD08, blue dots) and the GNIP database from both the Southern Hemisphere (black dots) and Northern Hemisphere (light grey dots). (b) Modeled relationship between δ18O of precipitation and surface temperature using our linear scaling, colored by initial evaporation air temperature (in C).


A3.3 Seasonality

Does seasonality in the hydrological cycle systematically bias climatological information in ice-core water-isotope records? While the difference between precipitation-weighted surface temperature and annual mean surface temperature is often discussed, this is not strictly the relevant comparison from the perspective of water-isotope ratios of precipitation. As discussed above, the critical comparison is between the annual mean surface temperature and the condensation-weighted temperature, integrated over both time and altitude. Our analysis of the MERRA-2 reanalysis data (Fig. A9) takes seasonal variation in precipitation and the vertical temperature profile into account. Differences in seasonality of condensation at different sites contributes to the spread around the central relationship. It is nevertheless useful to investigate potential bias in the precipitation-weighted surface temperature, since direct observations of Antarctic surface temperature are more common than full profiles of the atmospheric column.

The potential for precipitation weighting to bias annual average surface temperature depends on the phase angle between the seasonal cycles of precipitation and temperature. Only strong correlation or anticorrelation between the two cycles leads to persistent biasing. The potential for bias also depends on the ratio between stochastic and seasonal variability in both temperature and accumulation. If non-seasonal variance in accumulation is very large compared to the amplitude of the seasonal cycle in accumulation, for example, then the potential for bias is small. We examine seasonality in monthly surface temperature and snowfall over Antarctica in the ERA-Interim reanalysis as well as global precipitation in the MERRA-2 reanalysis. While the annual average Antarctic surface temperature and the precipitation-weighted surface temperature are often different in either reanalysis product, we find little systematic bias. Across the Antarctic the month-to-month and year-to-year variance in snowfall is large compared to the climatological seasonal cycle. The stochastic sampling of the seasonal cycle in surface temperature overwhelms the potential bias introduced by the average seasonality of precipitation. Further, the timing of the climatological annual cycle in snowfall varies across the continent, whereas the annual temperature cycle is quite coherent. The potential for seasonal bias thus varies dramatically between sites, even in the absence of the dominant stochastic sampling.

We compare the annual average surface temperature to the precipitation-weighted annual temperature at every grid point in the Southern Hemisphere. The mean bias for sites with typical Antarctic surface temperatures is less than 0.33 C, with the precipitation-weighted temperature being slightly colder on average. We find no systematic dependence of this bias on the surface temperature itself. While individual sites do show differences up to 4 C over the interval examined, our analysis does not suggest that such differences are persistent at a site. None of these analyses of monthly data take into account potential biases at the scale of individual precipitation events. The intermittency of Antarctic snowfall likely complicates the relationship between condensation-weighted and annual mean temperature at the seasonal and annual scale. At the same time, however, precipitation intermittency reduces potential seasonal biasing at climatological timescales by degrading any coherence in the seasonal cycles of accumulation and temperature.

Could seasonality in evaporation bias reconstructed source region T0? We examine the seasonality of Southern Hemisphere evaporation in the monthly MERRA-2 reanalysis by comparing the annual average over-sea surface temperatures to the evaporation-weighted annual temperatures. Between 35 and 65 S, the bulk of Antarctic moisture sources and evaporation-weighted surface temperatures are slightly colder than mean annual surface temperatures (the mean difference is 0.123 C, with 95 % of points between −0.25 and 0.5 C). South of the climatological sea ice zone, mean evaporation temperatures are a couple degrees warmer than mean annual surface temperatures on average, though our moisture tagging analysis suggests that these areas contribute at most a couple percent of the total moisture arriving at typical Antarctic sites.

Figure A11Partial derivatives of isotope state spaces with respect to surface temperature: (a) δ18OTs, (b) δDTs, (c) dxsTs, (d) dlnTs. Shading and contours in all panels is the slope in ‰ C−1.


A4 Tuning the Simple Water Isotope Model

We tune the Simple Water Isotope Model by adjusting the temperature dependence of the supersaturation of vapor over ice. Given insufficient observational and physical constraints, we parameterize the supersaturation as a linear function of temperature (Jouzel and Merlivat1984) as above, Si=a+b×T, set a=1, and tune the slope, b. The supersaturation has a strong influence on the kinetic fractionation (Eq. A12) and thus the relationship between δD and δ18O in vapor and precipitation. We tune the model to yield the observed relationship between δD and δ18O in global precipitation, rather than the relationship between δ values and environmental variables such as surface temperature.

Our target observational dataset includes water-isotope measurements of precipitation and surface snow from Antarctica and around the globe. The bulk of this compilation is that published by Masson-Delmotte et al. (2008). We include additional published surface snow and precipitation measurements from the GNIP database (IAEA2001), surface traverses at Dome A (Xiao et al.2013; Pang et al.2015), Dahe et al. (1994), and a 17Oxs compilation from Schoenemann et al. (2014). We also include previously unpublished measurements from a transect of snow pits and shallow firn cores across the main divide of the West Antarctic Ice Sheet. Samples from five sites were collected spanning 80 km across the ice flow divide in the 2012–2013 summer season. Samples were measured at IsoLab, University of Washington, Seattle, WA, USA. Measurement techniques are described in Markle et al. (2017). Measurements were made using laser spectroscopy (Picarro L2120-i analyzer). Data are reported relative to the VSMOW (Vienna Standard Mean Ocean Water) standard and normalized to SLAP.

The global relationship between δD and δ18O has an approximate slope of 8, as codified in the historical definition of the deuterium excess parameter (Dansgaard1964). However, the slope is not fundamental (Craig1961); as discussed in Sect. 1.2 the true observed relationship is nonlinear (Uemura et al.2012), as is the theoretical relationship even in the absence of kinetic fractionation (Markle et al.2017). Uemura et al. (2012) find an empirical fit between δD and δ18O in a global dataset for precipitation. They use a second-order polynomial fit, which is the basis for the logarithmic deuterium excess parameter (Uemura et al.2012; Markle et al.2017) (Eq. 4). From Eq. (A9), we can see that the theoretical relationship between δD and δ18O, given any amount of distillation, depends on the ratio of exp(Dαtot)exp(18αtot), where each αtot is itself a nonlinear function of temperature as outlined throughout Sect. A2.2. The ratio of exponential functions can be estimated to any arbitrary degree of accuracy with a polynomial function.

Our modern dataset includes several new sets of measurements in addition to those used in Uemura et al. (2012). We find similar coefficients in a second-order polynomial fit between δ18O and δD in our larger dataset compared to those found by Uemura et al. (2012): A=-29.2 and B=8.45 (see Eq. 4). Because these coefficients are not significantly different than those previously published and for consistency with that work, we use the coefficients found by Uemura et al. (2012) (A=-28.5 and B=8.47) to define a logarithmic deuterium excess parameter, dln. We find no benefit or justification for using higher-order fits to this dataset.

We run SWIM to produce an ensemble of temperature trajectories representing a wide range of possible evaporation and condensation temperatures (T0 varies from 0 to 28 C; Tc from 27 to −60C). We then compare the resulting cloud of modeled δ18O, δD, and dln, finding a second-order polynomial fit between the modeled δD and δ18O from the ensemble of temperature trajectories. We tune the model by iteratively adjusting the b value in the supersaturation parameterization to minimize the difference between the modeled and observed relationship between δD and δ18O (Eq. 4). This is easily visualized in a plot of δ18O and dln (e.g., Fig. A12a), as the average δ18O-to-dln relationship is flat in measured samples by definition.

Using the local closure assumption we find an optimal tuning of Si=1-b×T for b=0.00525C−1, as shown in Fig. A12. The observational data allow large ranges of the value of b to be rejected as shown in Fig. A13b and c; the fit coefficients of the resultant δD and δ18O are clearly irreconcilable with observations. While it is possible to optimize b as described above, there are limitations to this tuning procedure and the observational data may not allow discrimination within a small range of b values. In principle we should not expect a second-order polynomial fit between modeled δD and δ18O to be identical to the observed fit: the observational target data represent a variety of timescales from sub-seasonal to multiyear averages; the sites are neither evenly nor randomly distributed over the Antarctic continent, and the sites represented in the observational dataset have specific moisture source distributions, mean latitudes of evaporation, and evaporation temperatures which are not known a priori.

Because higher-elevation, colder Antarctic sites likely have more equatorward MSDs (Fig. 1), we should expect more depleted δ18O in the target data to be associated with slightly warmer T0 (that is, modeled results from a single value of T0 should not be strictly flat in the δ18Odln space). If we take the MSDs determined from the GCM experiments described earlier as representative and assuming the climatological meridional profile in surface air temperature, we should expect a 1–2 C increase in T0 between δ18O values of −40 ‰ and −55 ‰ in the observational dataset. This is a fairly small shift in mean T0 compared to the range of T0 that may contribute to a site but should give some downward curvature to model results of equal T0 at the coldest Tc values.

The appropriate weighting of model realizations with different T0 could vary within our tuning cost function, in principle, depending on the site of the target data. However, without knowing the true moisture source distribution and conditions for each sample in the target data a priori, assigning a single objective weighting scheme is difficult. We prefer values of b in which the more depleted observational data transect model realizations of slightly warmer T0, though this is not a strong constraint on the tuning. While one can reasonably reject most possible values of b (as in Fig. A13), we cannot justifiably discern between others within a small range (e.g., 0.005b0.0055). This gives a small but inherent uncertainty to the model tuning and in turn our temperature reconstructions. In spite of these limitations, the tuning procedure reproduces the observed isotope relationships well.

Figure A12Tuning SWIM with linear supersaturation parameterization. (a) The modeled δ18O and dln of precipitation (colored circles) for a range of condensation and evaporation temperatures. Color shading shows source region evaporation temperature in degrees Celsius (C). Black dots are the target dataset as described in the text. Modeled results are for the optimized supersaturation parameterization (Si=1-b×T, b=0.00525C−1) using the local closure assumption and the NCEP/NCAR reanalysis dataset for source region correlations. (b) Same as panel (a) for the modeled and target δ18O and δD of precipitation. (c) Same as panel (a) for the modeled and target δ18O and dxs of precipitation.


Figure A13Comparison of tuned and rejected supersaturation parameterizations, Si=1-b×T. (a) Tuned example: b=0.00525C−1. Same as in Fig. A12a. The modeled δ18O and dln of precipitation (colored circles) and target dataset (black dots). Color shading shows source region evaporation temperature in degrees Celsius (C). (b) Same as panel (a), but for a rejected tuning: b=0.003C−1. The modeled dln values curve upward strongly with δ18O, incongruently with the target data. (c) Same as panel (a), but for another rejected tuning: b=0.007C−1. The modeled dln values curve downward strongly with δ18O and incongruently with the target data.


The model tuning is not especially sensitive to which reanalysis dataset is used for correlations of the initial evaporation conditions or the season of evaporation. The model is, however, sensitive to which closure assumption is used.

While the linear, temperature-dependent parameterization of supersaturation is both simple and widely used, the physical processes determining supersaturation are complex. To understand the sensitivity of the model to this parameterization we also test a nonlinear parameterization of supersaturation, Si=a-b×T-c×T2. If the physical source of the high supersaturation at very cold conditions is related to the absence of condensation nuclei, the supersaturation may not linearly increase at very cold temperatures, as there could be diminishing returns as the atmosphere becomes increasingly clean. A small but positive c parameter that gradually decreases the slope in Si with decreasing temperature could be plausible. Only very small values of the second-order term, c, those of order 5×10-6C−2, are reconcilable with the observed δ18O-to-dln relationship. The modern data cannot readily distinguish whether the added complexity of the nonlinear parameterization is a better fit than the simple linear parameterization. For this reason, the linear parameterization is the most justifiable choice, though the uncertainty associated with this parameterization for temperature reconstructions will be examined below.

A5 Air parcel mixing within SWIM

We have so far assumed that the moisture-weighted average of a set of independent pseudo-adiabatic pathways, defined by a range of evaporation and condensation temperatures, can approximate the conditions experienced by precipitation falling at a polar site. We now aim to test the limits of this approximation and assess the influence of atmospheric mixing on the isotopic composition of air masses within the simple model framework.

The influence of air mass mixing during evaporation is considered in the discussion of the closure assumption above. Here we consider mixing during transport of air masses with different initial evaporation conditions, different condensation histories, and different temperature, moisture content, and isotopic values at the time of mixing. The central question is whether the processes of mixing can result in isotopic values of final precipitation that are significantly different than the moisture-weighted average of precipitation from unmixed pathways. There are three processes associated with mixing to consider.

  1. Non-uniqueness. Parcels that experienced different evaporation conditions can arrive at a condensation site of a given temperature with different isotopic values.

  2. Mixing-induced condensation. Mixing two saturated or undersaturated air parcels of different temperatures may result in an oversaturated mixed parcel due to nonlinearity in the Clausius–Clapeyron relationship. This process leads to additional condensation and fractionation (as well as warming due to latent heat release) and thus a more depleted isotopic value for a given temperature compared to the moisture-weighted average of unmixed pathways.

  3. Nonlinear mixing. For isobaric mixing of equal-massed air parcels, the final mixed temperature reflects their mass-weighted average (plus the effect of any latent heat release). However, the relative abundances of water isotopes mix with the moisture content of the air parcels rather than their total mass. Thus, while the temperatures have mixed linearly, the isotopic values of the resultant mixed parcel will be weighted nonlinearly toward those of the warmer and wetter parcel. The result is that the mixed parcel has a less depleted isotopic value for a given final temperature compared to the moisture-weighted average of the unmixed parcels.

Moisture-weighted differences between mixed and unmixed parcels only occur when air masses of different temperatures mix. Physically this may represent colliding fronts at synoptic scales. We consider two air parcels evaporated from identical starting conditions but which mix at different temperatures.

The influence of processes 2 and 3 is largest when the relative humidities of both parcels are at saturation, and the magnitude of the influences increases both as the difference in temperature between the two parcels increases and as the absolute temperature of the parcels increases. As temperature decreases the nonlinearity of the mixing of moisture approaches the linear (mass-weighted) mixing of temperature.

To assess the range of temperature differences associated with synoptic scales in the Southern Hemisphere, we examine the difference in daily mean 2 m air temperature from the ERA-Interim reanalysis (Dee et al.2011) over the Southern Hemisphere oceans. In summer, the day-to-day temperature differences have a standard deviation of less than 0.9 C, and in winter they have a standard deviation of less than 1.5 C (other reasonable metrics of synoptic-scale variability such as lagged 2 or 5 d temperature differences or grid-point-to grid-point differences are similar). While there is surely the potential for strong mixing for any given synoptic event, for the purposes of paleoclimate reconstruction we are interested in the long-term average of many storm events and thus the statistics of mixing generally.

We use the simple model to assess the range of final isotopic values of precipitation that can arise from two air parcels, which evaporated at the same initial conditions, then mixed at a range of different temperatures during transport. Air parcels are evaporated at specified initial air temperature (T0=10C), cooled, and randomly mixed at any combination of temperatures whose difference does not exceed a threshold, then cooled the remainder of the temperature pathway to −30C. We assume no preferential temperature of mixing and no preferential difference in temperature during mixing, though we assume a normally distributed probability of mixing with temperature differences up to 5 C (a high-end estimate based on the above analysis of daily temperature differences). This process is repeated 10 000 times to create a distribution of final δ18O (Fig. A14a), δD (Fig. A14b), and dln (Fig. A14c) values of precipitation at −30C, which is compared to the values from a parcel distilled along the same temperature pathway with no mixing (vertical red line in panels a–c, closed circle in panel d).

The resultant distributions are skewed and bimodal. The moisture-weighted means of the mixed distributions are shown as vertical black lines, while the unmixed final values are shown as vertical red lines. The means of the distribution are less depleted than the unmixed parcel owing largely to process 3. The influence of process 2 can also be seen in the additional peak at more depleted values. These distributions vary as a function of both T0 and Tc, though the differences in moisture-weighted means between mixed and unmixed parcels are relatively small and consistent.

The idealized tests above show the influence of mixing (at different temperatures) of air parcels that evaporate from identical source conditions. Perhaps more realistically, air parcels from different sources can mix at different air temperatures, in which case all three mixing processes above are important. This can act to broaden the distribution associated with a given moisture-weighted mean isotopic value. In Fig. A14e–h we show the distribution of the isotopic value of precipitation at −30C from a simulation of 10 000 randomly mixed air parcels as described above, except that the two parcels have two different initial evaporation temperatures of 5 and 15 C. We compare the distribution to the values associated with unmixed parcels originating at each evaporation temperature, as well as the moisture-weighted mean of the unmixed parcels. Differences between the mean isotopic values of mixed and unmixed parcels are less than 0.2 ‰ for δ18O, 1.5 ‰ for δD, and 0.01 ‰ for dln. Interestingly, although the distributions are broader in this scenario compared to the scenario in which the moisture parcels come from the same evaporative conditions, the differences between the moisture-weighted means of mixed and unmixed parcels are actually smaller. This is presumably because the skewed influence of process 3 (which drives the persistent bias above) contributes less to the total distribution. The relationship between δ18O and condensation temperature and the relationship between dln and evaporation temperature are similar whether or not mixing is present. Below we investigate the influence of mixing on our temperature reconstruction technique.

Figure A14Influence of air mass mixing during transport on water-isotope ratios of precipitation. (a–d) Air parcels from the same initial evaporation conditions (namely T0=10C) are distilled, stochastically mixed over a range different temperatures during transport, and distilled to a final precipitation temperature (−30C). Histograms of the final δ18O (a), δD (b), and dln (c) are shown, and the moisture-weighted means of those distributions are shown as the vertical black line. For comparison the final isotopic compositions from an unmixed pathway from 10 C to −30C are shown as the vertical red lines. The distribution in δ18Odln space of the 10 000 mixed pathways is shown in (d). (e–h) Same as for panels (a)(d) but for stochastically mixed air parcels arising from two different initial evaporation conditions (T0=5C and T0=15C). Moisture-weighted means of the mixed pathways are shown as vertical black lines, while the moisture-weighted means for equivalent unmixed pathways are in red. In (h) model results are colored by the moisture-weighted T0 resulting from the mixed pathways.


A6 Optimal coordinates for reconstruction technique

Consider a water sample with mean values of δ18O and δD and normally distributed uncertainties σ18 and σD. These uncertainties may arise from measurement uncertainty or in the mean δ value over some time or depth range represented by that sample. The dxs and dln values of the sample have corresponding σxs and σln, respectively, resulting from the propagation of σ18 and σD. How do these uncertainties influence the temperature reconstruction?

In Fig. A15 we examine the estimation of T0 using coordinates of δ18O and δD (Fig. A15a), δ18O and dxs (Fig. A15b), and δ18O and dln (Fig. A15c). The uncertainties in the position of the measurement along both the x axis (δ18O) and y axis (either δD, dxs, or dln) combine to give the total uncertainty in the position on the T0 surface, shown as the targets in Fig. A15a–c. The total combined uncertainty in the estimation of T0 is shown as probability density functions (PDFs) for each method in Fig. A15d. All estimates yield the same mean value of reconstructed T0; however, the widths of the probability density functions are different for each method. The δD method yields the broadest PDF and thus most uncertain reconstruction. While the PDFs for the dxs and dln reconstructions are similar, the dln reconstruction has a narrower PDF and thus more confident reconstruction. This is because the T0 isotherms are most separated along and most perpendicular to the dln axis. The advantage of the separation of isotherms along the dln axis is in part compensated for by the broadening of σln compared to σxs due to the propagation of uncertainties. However, the angle of the isotherms to the y axis is more important. Given a normal distribution of uncertainty along the y axis, perpendicular isotherms of the variable we wish to reconstruct will result in the narrowest possible distribution of that uncertainty across isotherms. If the angle of the isotherms deviates from perpendicular, as in the case with dxs at more depleted δ18O values, that uncertainty will be spread across a wider range of isotherms. The axis of the influence of T0 on dxs (and δD) is rotated with respect to its axis of variability in dxs.

Figure A15(a) Evaporation source temperature, T0, contours as a function of modeled δ18O and δD of precipitation. Uncertainty in δ18O and δD for an interval or sample is shown as PDFs of uncertainty along the respective axes. The intersection of these PDFs on the T0 surface results in a two-dimensional PDF in the reconstructed value of T0, shown as a target. (b, c) Same as panel (a) but for the evaporation source temperature projected onto the δ18O and dxs axes and the δ18O and dln axes, respectively. Uncertainties in δ18O and δD in panel (a) are propagated into the PDF on the dxs and dln axes in (b, c). (d) The uncertainty in the reconstructed evaporation source temperature (Tsource=T0) owing to the weighting of the combined two-dimensional PDFs from panels (a) (in blue), (b) (in red), and (c) (in purple).


This result is of course ultimately tied to the same reasons that dln provides a better qualitative proxy for source region changes than dxs. The initial imprint of the source conditions is better preserved in dln than dxs. The infidelity of the historical definition of the parameter is the result of nonlinear biases from the linear slope of the definition, the nonlinear nature of equilibrium fractionation, and the cumulative influence of kinetic fractionation during transport (Markle et al.2017).

A7 Reconstructions

In Fig. A16 we show the SWIM results (under our base assumptions) overlain with every pair of δ18O and dln measurements (corrected for changes in seawater δ18O) from the eight Antarctic deep ice-core records examined in this study.

Figure A16Inverted T0 and Tc surfaces as a function of modeled dln and δ18O of precipitation as in Fig. 7. Overlain are pairs of dln and δ18O measurements for eight different deep ice-core records, whose site locations are shown on the inset map. See Sect. 4.2 in the text for details on ice-core records.


We show the full reconstructions of source evaporation temperatures at all sites in Fig. A17, condensation temperatures in Fig. A18, and surface temperatures in Fig. A19. All sites have been resampled to even 200-year resolution for visual clarity. The light shading around the reconstructions represents the uncertainty in absolute temperature, while the darker shading represents the uncertainty in relative temperature variability as described in Sect. A9. Note that in some cases sample-to-sample variability is larger than the relative uncertainty, making the shading difficult to see at the scale plotted here.

Figure A17Source region evaporation air temperature reconstruction for all ice-core sites. The light shading around the reconstructions represents the uncertainty in absolute temperature, while the darker shading represents the uncertainty in relative temperature variability as described in Sect. A9.


Figure A18Condensation air temperature reconstruction for all ice-core sites. The light shading around the reconstructions represents the uncertainty in absolute temperature, while the darker shading represents the uncertainty in relative temperature variability as described in Sect. A9.


Figure A19Surface air temperature reconstruction for all ice-core sites. The light shading around the reconstructions represents the uncertainty in absolute temperature, while the darker shading represents the uncertainty in relative temperature variability as described in Sect. A9.


A8 Correlation of nonlinear and linear reconstruction techniques

In Fig. A20 we show the difference between the nonlinear reconstructions for all core sites and reconstructions based on linearization around Holocene and glacial conditions as described in the text (Sect. 4.3.1). We also show the difference in the reconstructions using the Holocene and glacial linearizations.

Figure A20Difference in site surface temperature and evaporation source air temperature using SWIM results and different reconstruction techniques. (a, b) Difference between the nonlinear reconstruction technique and linearization of SWIM results around Holocene conditions. (c, d) Difference between the nonlinear reconstruction technique and linearization of SWIM results around glacial conditions. (e, f) Difference in reconstructions using linearizations of SWIM results around both Holocene and glacial conditions.


We compare the nonlinear temperature reconstructions of all eight ice-core sites to linearized reconstructions using SWIM results for Holocene conditions as described in Sect. 4.3.1. We interpolate all records to even time spacing and compute correlation matrices amongst cores for T0 and Ts using the linear and nonlinear technique (Fig. A21). Reconstructed site surface temperatures, Ts, are extremely well correlated amongst all cores using either technique, though there is marginal improvement in correlation using the nonlinear technique. In the case of evaporation temperatures, T0, there is dramatic improvement in coherence amongst the records when using the nonlinear technique. The increase in shared variance (R2) explained using the nonlinear technique is shown in Fig. A22. Note that the largest increase in shared variance is associated with the Siple record. This makes sense given the conditions of that site compared to the others and the patterns of partial slopes in Figs. 4 and 5.

By accounting for the fundamental nonlinearities in water-isotope distillation we are able to reveal more coherent underlying climate signals in source region temperatures, which are otherwise obscured by linear temperature reconstruction techniques. For analogous reasons, Markle et al. (2017) argued that the logarithmic deuterium excess parameter dln is a more faithful qualitative proxy for source region conditions than the linearly define dxs. Compare the correlation matrices of the excess parameters in Fig. A23. The nonlinearly reconstructed T0 and dln parameters share the same correlation pattern amongst the ice cores and show substantially more coherence than either linearly reconstructed T0 or dxs. The correlation pattern of dln and dxs between all core sites (Fig. A24) reveals how nonlinear effects alter the traditionally defined dxs at the coldest Antarctic temperatures. The broad change from positive to negative correlation of dln to dxs across sites is a reflection of the change in sign of dxsTc as a function of Tc.

Figure A21Correlation (R) matrices for temperature reconstructions of all core sites. (a) Correlation among reconstructions of T0 using a linearization of SWIM results for isotopic conditions of the Holocene. (b) Correlation among reconstructions of T0 using the full nonlinear SWIM results. (c) Correlation among reconstructions of Ts using a linearization of SWIM results for isotopic conditions of the Holocene. (d) Correlation among reconstructions of Ts using the full nonlinear SWIM results. All records are ordered by their approximate modern surface temperature.


Figure A22The increase in R2 (shared variance) of the nonlinear reconstructions of T0 for all sites over the linear reconstructions (linearized around Holocene conditions). All records are ordered by their approximate modern surface temperature.


Figure A23Correlation matrices for (a) δ18O, (b) δD, (c) dxs, and (d) dln among all core sites. All records are ordered by their approximate modern surface temperature.


Figure A24Correlation matrix between dxs and dln between all core sites. All records are ordered by their approximate modern surface temperature.


A9 Temperature reconstruction uncertainty

In this section we investigate uncertainty in our temperature reconstructions by examining the sensitivity of our results to assumptions and parameterizations in the model. We can compare reconstructed Tc and T0 from a set of δ18O and dln measurements using multiple iterations of the model in which the value of a parameter or an underlying assumption has been varied. In Fig. A25 we show Tc and T0 reconstructions for the WDC record (Markle et al.2017) arising from a number of model parameters and assumptions, which are discussed below. Because our reconstruction technique takes into account nonlinearities, differences in reconstructed temperatures may have mean offsets and may have differences in variability that vary as a function of δ18O and dln. Thus, uncertainty arising from a given parameter may vary between ice-core sites. In general, varying a parameter in the model results in patterns of the partial slopes in δ18O and dln with Tc and T0 that are similar, but shifted in the Tc and T0 space. A consequence of this is that the uncertainty in absolute values of reconstructed Tc and T0 is generally larger than uncertainty in their relative variability.

Figure A25Temperature reconstruction sensitivity to model parameterizations. Tc and T0 reconstructions using WDC δ18O and dln (resampled to 10-year resolution) as well as SWIM runs. Black lines show the base model in all panels, while colored lines are results in which model parameters and assumptions are varied: (a–b) evaporation condition correlations based on NCEP (black) and ERA (red) reanalysis data; (c–d) global (blue) and local (black) closure assumption during evaporation. (e–f) A range of values for 18αdiff; (g–h) a range of values for the b parameter in the supersaturation parameterization, as well as a nonlinear parameterization (c) as described in the text; (i–j) several versions of the precipitation parameterization in which all moisture is removed above saturation (sat), all moisture is removed above initial RH (RH = RH0, constant RH along path), and all moisture is removed above fixed 80 % or 90 % RH.


It is useful to distinguish between uncertainty in the true value of a parameter in the modern climate and the possibility that the effective value may change as a function of climate. Further, not all sources of uncertainty are independent. Varying the value of some parameters may require retuning the model before calculating the isotope state spaces. By ignoring this we risk conflating uncertainty in the reconstruction with bias in the reproduction of the modern mean state.

A9.1 Sensitivity to model parameters

There is uncertainty in our reconstructions associated with the tuning procedure. While we can constrain the possible values of the b parameter in the supersaturation function by comparison to modern data, variations within a small range should not be ruled out given the imperfect constraint of modern observations. In Fig. A25 we show the resultant uncertainty in the temperature reconstruction from uncertainty in the supersaturation parameterization (b=0.0051 to 0.0054C−1). Uncertainty arising from other aspects of the distillation scheme, such as the value of the diffusive fractionation factors during transport, is encapsulated by the tuning uncertainty since adjusting those parameters requires retuning the model.

Aspects of the initial evaporation scheme introduce uncertainty into our reconstructions. The value of αdiff18 during evaporation is important in setting the initial isotopic values of vapor. While we find the value 1.009 to give the best fit to modern observations, values within a small range may be defensible (Fig. A6). The local closure assumption used in the evaporation scheme has known limitations (Risi et al.2010), representing an end-member scenario for possible evaporative conditions. While less applicable to past climate mean states, the global closure assumption provides an extreme test of the model's sensitivity. Using the global rather than local closure assumption can lead to differences in reconstructed absolute T0 up to 1.5 C for the WDC record, while differences in absolute Tc are smaller (≤1C). Relative variability in T0 and Tc is similar when using either closure assumption and ≤0.3C.

We also examine the influence of the source relative humidity parameterization on our temperature reconstructions. In our base model we use climatological correlations to determine an initial relative humidity given an initial air temperature; colder surface air temperatures over the ocean are associated with slightly higher relative humidity. We show the difference in reconstructions due to using either the NCEP or ERA-Interim reanalysis. We can also ask how different our reconstructed Tc and T0 in WDC would be if we used fixed mean values of initial relative humidity rather than values that depend on T0. These differences are not true uncertainties in the reconstruction as variable surface relative humidity is a more physically defensible choice than a fixed relative humidity, though these tests serve to demonstrate the robustness of the reconstruction to model assumptions. In a similar vein we can examine the sensitivity of the model to our precipitation parameterization and the potential choices of that parameterization discussed in Sect. A1.2.

A9.2 Influence of mixing on temperature reconstruction

We assess the potential influence of atmospheric mixing on our temperature reconstruction framework by comparing the maps of T0 and Tc as functions of δ18O and dln to maps produced by a large ensemble of model runs that incorporate stochastic mixing. We consider a range of final condensation temperatures from −20 to −50C. We generate random pairs of pseudo-adiabatic cooling pathways ending at every value of Tc and random values of T0 pulled from a normal distribution (with a mean of 12 C and standard deviation of 4 C) similar to the modern Antarctic moisture source distributions. Air parcels cooled down these two paths are stochastically mixed at points along the path and cooled to the final Tc. To mix, parcel temperatures must be above an absolute threshold temperature (−15C) and have a relative difference within 5 C, as described above. This results in a conservative estimate of the influence of mixing: mixing at lower temperatures reduces the average difference between mixed and unmixed pathways since the effects of mixing are larger when absolute humidity is higher. We take 50 random mixtures from each of 2×50 random cooling pathways for each value of Tc between −20 and −50C in increments of 0.1 C (a total of 1.5×104 unmixed and 7.5×105 mixed cooling paths). We then interpolate the results of both the mixed pathways and the moisture-weighted averages of the unmixed pathways to create maps of T0 and Tc as functions of δ18O and dln (at a resolution of 0.1 ‰, Fig. A26). Due to the stochastic mixing these maps are unevenly populated. The potential influence of mixing on our reconstruction technique can be seen in the difference between the mean Tc and T0 maps resulting from the mixed (denoted with the subscript M) and unmixed pathways (Fig. A26). In Fig. A26 we show the histograms of ΔTc=Tc,M-Tc and ΔT0=T0,M-T0. We test a range of mixing and threshold values in multiple Monte Carlo simulations. In all cases the mean values of Tc,M-Tc and T0,M-T0 are very near zero (<|±0.06C| for Tc, and <|±0.02C| for T0). The spread of the histograms in Fig. A26 represents the inherent uncertainty in our reconstruction technique when mixing is neglected. This uncertainty is less than ±0.2C for Tc and ±0.35C for T0 in all tests. We find similar results when using isobaric rather than pseudo-adiabatic cooling pathways. Including moisture sources with T0<0C in our mixing analysis has no significant influence on the mean difference between the weighted mean maps of either Tc or T0, though expanding the range of moisture sources to include T0<0C does increase the range of T0,M-T0 by over a degree, in agreement with the analysis of unmixed pathways.

Figure A26The influence of air parcel mixing on modeled isotope state spaces. (a, e) The water-isotope ratios of precipitation resulting from 7.5×105 stochastically mixed distillation pathways, colored by Tc and T0, respectively. (b, f) The water-isotope ratios of precipitation resulting from the corresponding 1.5×104 unmixed distillation pathways, colored by Tc and T0, respectively. (c, g) The difference in Tc and T0, respectively, between the mixed and unmixed pathways as a function of the δ18O and dln space. (d, h) Histograms of the differences, ΔTc and ΔT0, from all points in the δ18O and dln space, resulting from mixing.


A9.3 Combined uncertainty estimates

To calculate the total uncertainty in our temperature reconstructions, we examine the combined influence of the major independent sources of uncertainty. These include tuning via the supersaturation function, the closure assumption, mixing in the atmosphere during transport, the precipitation parameterization, the diffusive fractionation factor during evaporation, and the relationship between initial air temperature and relative humidity. We calculate the absolute uncertainty for each Tc and T0 reconstruction interpolated at each pair of δ18O and dln measurements as the absolute difference in reconstructions arising from perturbations to parameter values or assumptions. To estimate the uncertainty in relative temperature changes we subtract the mean of each reconstruction before calculating differences due to parameter perturbations.

We estimate the uncertainty due to model tuning as the mean absolute difference from the base scenario for reconstructions using values of b=0.0051 to 0.0054 in the supersaturation function. Likewise the impact of uncertainty in the value of the diffusive fractionation factor is estimated as the mean absolute difference of reconstructions using 18adiff=1.009±0.001. We estimate the uncertainty introduced by the precipitation parameterization as the mean absolute difference from the base scenario of reconstructions using each of the alternate assumptions outlined in Sect. A1.2, applied symmetrically to the base scenario. We estimate the uncertainty arising from our assumed relationship between T0 and RH0 as the mean absolute difference in reconstructions using climatological fits from the NCEP/NCAR and ERA-Interim reanalysis.

Based on the tests in Sect. A2.1, a conservative estimate of the uncertainty arising from mixing at the evaporation source is half the absolute mean difference in reconstructions employing the local and global closure assumptions, applied symmetrically about the base scenario (local closure). Because of the stochastic nature of our atmospheric mixing simulations, our estimates of the Tc and T0 differences are nonuniform and unevenly populated, making interpolation in the δ18Odln space challenging (see Sect. A5). We thus take a conservative estimate of the absolute and relative uncertainty introduced by mixing during transport as the mean and standard deviation of the differences in the mixed and unmixed reconstructions across the entire state space, respectively (see Fig. A14).

We add the uncertainty from each independent source in quadrature as functions of δ18O and dln symmetrically around the base scenario. Finally, we include the additional uncertainty in our estimates of relative Ts variability arising from the Tc-to-Ts relationship outlined in Sect. A3.2. We use the mean absolute difference of reconstructions using Tc0.69±0.02CCTs to estimate this uncertainty, which is added in quadrature to the above uncertainties in Ts. An example of this spread of uncertainty in the WDC Ts reconstruction is shown in Fig. A27. Reconstructions of Ts, Tc, and T0 for several major ice-core records along with the combined relative uncertainty are those reconstructions is shown in Fig. 9.

Figure A27Sensitivity of the Ts reconstruction for WDC to the relationship between Tc and Ts. Our base scenario is shown in black (Tc0.69CCTs), while the spread associated with a ±0.02CC uncertainty in the scaling factor is shown in red and gold. For comparison the results from a weaker slope of 0.65CC are shown in blue. This range of scaling factors has little impact on the reconstructed temperature history and may be difficult to see.


Figure A28Annual mean temperature of moisture source regions from moisture-tagged GCM experiments. (a) Annual mean zonal-mean surface air temperature vs. latitude (blue) overlain with annual mean MSDs for several ice-core sites, colored by MSD-weighted annual mean air temperature. (b) MSD-weighted annual mean air temperature for every model grid point over the Antarctic. Note that this is not strictly the evaporation-weighted air temperature; it is weighted by location but not necessarily by time and may thus be biased to be too low.

Figure A29Isotope model results colored by initial evaporation air temperature T0 (C) using base model assumptions. For initial evaporation air temperatures below 0 C there are non-unique results in the δ18O and dln space.


Figure A30Temperature reconstructions for WDC accounting for below-zero evaporation air temperatures. (a) Reconstructions of Tc under base conditions (blue; no contribution from T0<0C), with the contribution from T0<0C weighted by the final rs (red) and with a 90 % contribution from T0>0C as well as 10 % from T0<0C (gold). (b) Same as in (a) but for reconstructions of T0.


A9.4 Uniqueness and source temperature

All Antarctic sites have mean initial evaporation air temperatures above 0 C, according to the moisture source distributions from water-tagged GCM experiments (Fig. A28). In fact, 85 %–95 % of all moisture that arrives at Antarctic sites in our modeling initially evaporates from locations with annual average surface air temperatures above freezing. The relatively small but nonzero contribution of moisture from evaporation temperatures below freezing poses an interesting challenge to our temperature reconstruction method. While it is widely known that there is not a unique value of δ18O for every value of condensation temperature owing to the influence of evaporation temperature, there are not necessarily unique pairs of δ18O and dln for every pair of Tc and T0 if T0 can be both above and below freezing. The Tc and T0 surfaces fold over on themselves in the δ18O and dln space for values of T0 below 0 C. An example of such a folded surface is shown in Fig. A29. Given a lack of isotopic vapor measurements for evaporation air temperatures much below 0 C and that our evaporation scheme is not well calibrated for such conditions, these results are purely illustrative.

Those caveats aside, we investigate the sensitivity of our temperature reconstructions to this potential non-uniqueness in the water-isotope state spaces. We run SWIM through a large field of T0 values from −28 to 28 C. We can in principle resolve the non-uniqueness problem by combining reconstructions from either side of the folded surface based on the contribution of total moisture represented by each pair of non-unique paths. Knowing that below-zero moisture sources contribute far less to the total moisture reaching Antarctic sites (Fig. A28), we examine two reasonable methods of moisture-weighting the reconstructions. In the first we simply weight each pair of reconstructed temperatures by the final mixing ratio (rs(eff)) of each modeled distillation path. This approach has the advantage of allowing contributions to vary with temperature and thus mean climate and leads to roughly 10 %–20 % contributions from below 0 C moisture sources to modern Antarctic sites (using a local closure assumption). However, this approach ignores the influence of dynamics and topographic–energetic isolation in determining Antarctic moisture source distributions, ultimately overestimating contributions from below 0 C moisture sources to higher-elevation, colder sites compared to our GCM-based MSD estimates. In this rs(eff)-weighted scheme, higher Antarctic sites have a relatively greater contribution from colder moisture sources than warmer sites owing to the curvature in the Clausius–Clapeyron relationship. Our moisture tagging analysis and previous studies (e.g., Bailey et al.2019) suggest, however, that transport dynamics should lead to the opposite relationship. An alternate approach is to specify fixed contributions from above and below 0 C moisture sources (e.g., 90 % T0>0C, 10 % T0<0C). While these average relative contributions are based on our moisture tagging analysis, we do not specify contributions as a function of site elevation or mean climate. Reconstructions based on these approaches (both calculated and specified moisture weighting) are shown in Fig. A30 for the WDC record. Considering the non-uniqueness leads to very small differences in reconstructed Tc and T0 variability: the standard deviation of differences in reconstructed Tc is less than 0.07 C using either method and less than 0.19 C for T0. Attempting to account for this non-uniqueness does, however, lead to persistent mean offsets in absolute temperature; in particular, we find colder absolute values of reconstructed T0 for all ice-core sites.

While these results are interesting, this attempt to account for non-uniqueness likely does not actually improve the absolute temperature reconstructions. Given the shape of the folded temperature surfaces in the modeled δ18O and dln space, as well as the actual values of δ18O and dln in ice-core measurements, the model must extrapolate to extremely cold T0 values for the below 0 C side of the folded surface. These values of T0 are far colder (>10C colder) than realistic evaporation temperatures likely to contribute moisture to high Antarctic sites given energetic constraints (Bailey et al.2019) and our moisture tagging GCM experiments. Further, as stated above, our evaporation scheme is not well calibrated to such evaporation conditions. Near-surface relative humidity in particular is not well constrained by our climatological correlations in these circumstances. Our model likely underestimates the depletion of the initial evaporate in these circumstances, meaning that the reconstruction solves for a very cold T0 when a much warmer one (and perhaps a reduced RH0) is actually correct. The net result of these considerations is that the analysis above should represent a quite conservative estimate of the influence of non-uniqueness on our temperature reconstructions and their relative variability; the real effect is likely much smaller though difficult to quantify precisely.

A10 Comparison to previous reconstructions

We next reconstruct site and source temperatures for four East Antarctic ice-core records and compare to previously published linear reconstructions. We use records of δ18O and δD measurements (and calculate dxs and dln) from the Vostok (Jouzel et al.1997; Uemura et al.2012), EPICA Dome C (EDC) (Stenni et al.2004, 2010), EPICA Dronning Maud Land (EDML) (Stenni et al.2010), and Dome Fuji records (Uemura et al.2012). After seawater correction, we use the ice-core δD and dxs for the linear reconstruction and δ18O and dln for the nonlinear reconstruction. The linear reconstruction parameters from several studies are compiled by Uemura et al. (2012) (see Tables 1 and 2 in Uemura et al.2012). Previous reconstruction techniques solve for the source temperature, Tsource, equivalent to our evaporation temperature, T0, and for the site surface temperature, Tsite. We convert our reconstructed condensation temperatures, Tc, to surface temperatures following the method in Sect. A3.2.

A comparison of relative changes in site and source temperatures is shown in Fig. 11. The nonlinear reconstruction results of this study are shown in black, while published linear inversions for each core are shown in color. The differences between the results of this study and the previous temperature reconstructions arise from differences between the linear and nonlinear reconstruction techniques as well as differences in the underlying water-isotope models used for the estimation of scaling relationships. In many cases, the previously published linear inversions overestimate changes in both site and source temperature compared to the nonlinear reconstruction.

The overestimation of reconstructed temperature change by the linear reconstruction makes physical sense. The largest source of nonlinearities in the water-isotope to temperature relationships is the deuterium excess parameters, dxsTc and dxsT0. If one assumes these slopes are linear over a given range in T0 and Tc when in reality they are nonlinear, one will attribute a given change in Δdxs to a larger change in temperature than is actually required. This overestimate of the required temperature change will be distributed across the reconstructed site and source temperatures in proportion to the values of the β and γ parameters. The same reasoning is true for nonlinearities in the relationships between δD or δ18O and the temperature boundary conditions, though the nonlinearities in these slopes are much smaller.

The residuals between relative temperature change in the nonlinear and previous linear reconstructions are shown in Fig. 12. Residuals in the site temperature reconstructions are on the order of ±2C (Fig. 12a). The residuals are not random but rather correlated with the reconstructions themselves, pointing to nonlinear biases.

The previous reconstructions use a different scaling between the surface and condensation than that used in this study (see Sect. A3.2). However, the differences between the nonlinear reconstruction and the linear reconstructions do not arise solely because of this different surface–condensation temperature scaling. The residuals between reconstructed condensation temperatures are shown in Fig. 12b. These differences are somewhat damped compared to those of the surface temperatures owing to different assumed slopes in the condensation to surface temperature relationship but are of similar magnitude, and the time series of the residuals are again correlated with the reconstructions themselves.

The residuals between the reconstructed evaporation temperature anomalies (Fig. 12c) have a large spread ranging from about +3 to −5C. While the magnitudes of source temperature residuals are comparable to those of site temperature, they are far more significant, representing from 50 % to over 200 % of the total reconstructed variability in the source temperature.

The residuals between the reconstructed evaporation temperature anomalies (Fig. 12c) have a large spread ranging from about +3 to −5C. As discussed above the largest source of potential biases is the deuterium excess relationships to temperature and should be greatest in the reconstruction of source temperatures. While the magnitude of source temperature residuals is comparable to those of site temperature, they are far more significant, representing between 50 % and over 200 % of the total reconstructed variability in the source temperature. This is related to the issues surrounding the qualitative interpretation of source region changes from dxs versus dln (Markle et al.2017; Uemura et al.2012) (see Sect. 1.2) and ultimately a consequence of the same distillation effects.

A11 Three-parameter reconstructions

In the approach outlined above, we have considered the boundary conditions Tc and T0 to be the only independent input variables. In particular, we have assumed that the source region relative humidity, RH0, is a dependent variable whose value is not fixed but determined by climatological correlation with T0. Most previous linear reconstructions have calculated scaling factors based on fixed values of RH or the average variation in RH over some range (Uemura et al.2012; Winkler et al.2012).

We can relax the assumption that RH0 is dependent on T0 and reconstruct three independent climate variables (Tc, T0, and RH0) if we have three independent constraints. While δ18O and dln alone are not sufficient, 17Oxs=δ17O-0.528δ18O (Landais et al.2008) can in principle provide the necessary additional information. We can allow Tc, T0, and RH0, to all vary as independent variables, defining a three-dimensional parameter space through which SWIM is run to produce three-dimensional isotope state spaces.

While promising, this method currently has practical limitations. Our model does not reproduce the observed 17Oxs-to-δ18O relationship in Antarctic precipitation to sufficient precision to offer useful constraints. This may be a consequence of model limitations such as missing physical processes. Alternatively (or additionally) uncertainties in the absolute values of 17Oxs in Antarctic precipitation may be too large to offer useful discrimination amongst variations in Tc, T0, and RH0 (Schoenemann et al.2014).

The 17Oxs of Antarctic precipitation in our model is sensitive to the supersaturation, diffusivities, and other parameters driving kinetic fractionation. Both small changes in the supersaturation parameterization and uncertainties in the absolute value of 17Oxs lead to large changes in the absolute value of reconstructed source region conditions (T0 and RH0). It is worth noting that absolute values of 17Oxs are 3 orders of magnitude smaller than the deuterium excess. Further, preliminary testing suggests that there may be significant non-uniqueness to address; that is, a position in the three-parameter space, defined by δ18O, dln, and 17Oxs, does not necessarily lead to unique values of the boundary conditions.

This general approach is scalable. Additional quantities that are both influenced by the environmental pathway and measurable in an ice core, for example accumulation rate (Fudge et al.2016), water-isotope diffusion lengths (Johnsen et al.2000), or the concentration of aerosols (Markle et al.2018), may be added to the model. These additional proxies can allow for the reconstruction of additional independent variables and the relaxation of assumptions. Alternatively it may be possible to use the same approach to optimize model parameters like the supersaturation. We leave this task to future work.

Code and data availability

The temperature reconstructions, underlying data, and Simple Water Isotope Model code are available through a public GitHub repository (doi:, Markle and Steig2022;, Markle2022) and the United States Antarctic Research Program (USAP) Data Center. The ice-core data are all already available through links provided in the primary papers cited for each dataset.

Author contributions

BRM conceived of the study and, with EJS, designed the model and analyses. BRM and EJS wrote the paper.

Competing interests

The contact author has declared that neither they nor their co-author has any competing interests.


Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.


Bradley R. Markle was supported in part by the Stanback Postdoctoral Fellowship. The authors wish to thank Hans Christian Steen-Larsen, Spruce Schoenemann, Emma Kahle, Vasileios Gkinis, Marina Dütsch, and Peter Blossey for helpful discussions. We would like to thank Robert Hellmann and two anonymous reviewers.

Financial support

This research has been supported by the National Science Foundation (grant nos. 1043092, 1141839, and 1443105).

Review statement

This paper was edited by Carlo Barbante and reviewed by two anonymous referees.


Bailey, A., Singh, H. K., and Nusbaumer, J.: Evaluating a Moist Isentropic Framework for Poleward Moisture Transport: Implications for Water Isotopes over Antarctica, Geophys. Res. Lett., 46, 7819–7827,, 2019. a, b, c, d

Bakhshaii, A. and Stull, R.: Saturated Pseudoadiabats-A Noniterative Approximation, J. Appl. Meteorol. Climatol., 52, 5–15, 2013. a

Barkan, E. and Luz, B.: High precision measurements of 17O/16O and 18O/16O ratios in H2O, Rapid Commun. Mass Spect., 19, 3737–3742, 2005. a

Barkan, E. and Luz, B.: Diffusivity fractionations of H216O/H 217O and H216O/H218O in air and their implications for isotope hydrology, Rapid Commun. Mass Spect., 21, 2999–3005, 2007. a, b

Benetti, M., Steen-Larsen, H. C., Reverdin, G., Sveinbjörnsdóttir, Á. E., Aloisi, G., Berkelhammer, M. B., Bourlès, B., Bourras, D., De Coetlogon, G., Cosgrove, A., Faber, A.-K., Grelet, J., Hansen, S. B., Johnson, R., Legoff, H., Martin, N., Peters, A. J., Popp, T. J., Reynaud, T., and Winther, M.: Stable isotopes in the atmospheric marine boundary layer water vapour over the Atlantic Ocean, 2012–2015, Sci. Data, 4, 1–17, 2017. a, b, c

Bintanja, R. and Van de Wal, R.: North American ice-sheet dynamics and the onset of 100,000-year glacial cycles, Nature, 454, 869–872, 2008. a

Bosilovich, M. G., Robertson, F. R., Takacs, L., Molod, A., and Mocko, D.: Atmospheric water balance and variability in the MERRA-2 reanalysis, J. Climate, 30, 1177–1196, 2017. a

Braconnot, P., Harrison, S. P., Kageyama, M., Bartlein, P. J., Masson-Delmotte, V., Abe-Ouchi, A., Otto-Bliesner, B., and Zhao, Y.: Evaluation of climate models using palaeoclimatic data, Nat. Clim. Change, 2, 417–424, 2012. a

Brook, E. J., White, J. W., Schilla, A. S., Bender, M. L., Barnett, B., Severinghaus, J. P., Taylor, K. C., Alley, R. B., and Steig, E. J.: Timing of millennial-scale climate change at Siple Dome, West Antarctica, during the last glacial period Siple Dome, West Antarctica, during the last glacial period, Quaternary Sci. Rev., 24, 1333–1343, 2005. a

Buizert, C., Sigl, M., Severi, M., Markle, B. R., Wettstein, J. J., McConnell, J. R., Pedro, J. B., Sodemann, H., Goto-Azuma, K., Kawamura, K., Fujita, S., Motoyama, H., Hirabayashi, M., Uemura, R., Stenni, B., Parrenin, F., He, F., Fudge, T. J., and Steig, E.: Abrupt ice-age shifts in southern westerly winds and Antarctic climate forced from the north, Nature, 563, 681–685, 2018. a, b, c, d

Buizert, C., Fudge, T., Roberts, W. H., Steig, E. J., Sherriff-Tadano, S., Ritz, C., Lefebvre, E., Edwards, J., Kawamura, K., Oyabu, I., Motoyama, H., Kahle, E. C., Jones, T. R., Abe-Ouchi, A., Obase, T., Martin, C., Corr, H., Severinghaus, J. P., Beaudette, R., Epifanio, J. A., Brook, E. J., Martin, K., Chappellaz, J., Aoki, S., Nakazawa, T., Sowers, T. A., Alley, R. B., Ahn, J., Sigl, M., Severi, M., Dunbar, N. W., Svensson, A., Fegyveresi, J. M., He, C., Liu, Z., Zhu, J., Otto-Bliesner, B. L., Lipenkov, V. Y., Kageyama, M., and Schwander, J.​​​​​​​: Antarctic surface temperature and elevation during the Last Glacial Maximum, Science, 372, 1097–1101, 2021. a, b

Ciais, P. and Jouzel, J.: Deuterium and oxygen 18 in precipitation: Isotopic model, including mixed cloud processes, J. Geophys. Res.-Atmos., 99, 16793–16803, 1994. a, b, c

Connolley, W.: The Antarctic temperature inversion, Int. J. Clim., 16, 1333–1342, 1996. a, b, c, d

Craig, H.: Isotopic variations in meteoric waters, Science, 133, 1702–1703, 1961. a, b, c

Craig, H. and Gordon, L. I.: Deuterium and oxygen 18 variations in the ocean and the marine atmosphere, in: Stable Isotopes in Oceanographic Studies and Paleotemperatures, Conferences in Nuclear Geology, 26–30 July 1965, Spoleto, Italy, edited by: Tongiorgi, E., 1965. a

Criss, R. E.: Principles of stable isotope distribution, Oxford University Press, ISBN 9780195117752, 1999. a, b, c, d, e, f, g

Cuffey, K. M., Clow, G. D., Steig, E. J., Buizert, C., Fudge, T., Koutnik, M., Waddington, E. D., Alley, R. B., and Severinghaus, J. P.: Deglacial temperature history of West Antarctica, P. Natl. Acad. Sci. USA, 113, 14249–14254, 2016. a, b, c

Dahe, Q., Petit, J., Jouzel, J., and Stievenard, M.: Distribution of stable isotopes in surface snow along the route of the 1990 International Trans-Antarctica Expedition, J. Glaciol., 40, 107–118, 1994. a

Dai, A.: Recent climatology, variability, and trends in global surface humidity, J. Climate, 19, 3589–3606, 2006. a

Dansgaard, W.: Stable isotopes in precipitation, Tellus, 16, 436–468,, 1964. a, b, c, d, e, f, g, h

Dansgaard, W., Johnsen, S. J., Møller, J., and Langway, C. C.: One thousand centuries of climatic record from Camp Century on the Greenland ice sheet, Science, 166, 377–380, 1969. a

Dee, D., Uppala, S., Simmons, A., Berrisford, P., Poli, P., Kobayashi, S., Andrae, U., Balmaseda, M., Balsamo, G., Bauer, P., Bechtold, P., Beljaars, A. C. M., van de Berg, L., Bidlot, J., Bormann, N., Delsol, C., Dragani, R., Fuentes, M., Geer, A. J., Haimberger, L., Healy, S. B., Hersbach, H., Hólm, E. V., Isaksen, L., Kållberg, P., Köhler, M., Matricardi, M., McNally, A. P., Monge-Sanz, B. M., Morcrette, J.-J., Park, B.-K., Peubey, C., de Rosnay, P., Tavolato, C., Thépaut, J.-N., and Vitart, F.: The ERA-Interim reanalysis: Configuration and performance of the data assimilation system, Q. J. Roy. Meteorol. Soc., 137, 553–597, 2011. a, b

EPICA Community Members: One-to-one coupling of glacial climate variability in Greenland and Antarctica, Nature, 444, 195–198,, 2006. a

Epstein, S., Sharp, R. P., and Goddard, I.: Oxygen-isotope ratios in Antarctic snow, firn, and ice, The J. Geol., 71, 698–720, 1963. a

Fudge, T., Markle, B. R., Cuffey, K. M., Buizert, C., Taylor, K. C., Steig, E. J., Waddington, E. D., Conway, H., and Koutnik, M.: Variable relationship between accumulation and temperature in West Antarctica for the past 31,000 years, Geophys. Res. Lett., 43, 3795–3803, 2016. a

Gelaro, R., McCarty, W., Suárez, M. J., Todling, R., Molod, A., Takacs, L., Randles, C. A., Darmenov, A., Bosilovich, M. G., Reichle, R., Wargan, K., Coy, L., Cullather, R., Draper, C., Akella, S., Buchard, V., Conaty, A., da Silva, A. M., Gu, W., Kim, G.-K., Koster, R., Lucchesi, R., Merkova, D., Nielsen, J. E., Partyka, G., Pawson, S., Putman, W., Rienecker, M., Schubert, S. D., Sienkiewicz, M., and Zhao, B.​​​​​​​: The modern-era retrospective analysis for research and applications, version 2 (MERRA-2), J. Climate, 30, 5419–5454, 2017. a, b

Gkinis, V., Holme, C., Kahle, E. C., Stevens, M. C., Steig, E. J., and Vinther, B. M.: Numerical experiments on firn isotope diffusion with the Community Firn Model, J. Glaciol., 67, 450–472, 2021. a

Gonfiantini, R.: Oxygen isotope variations in Antarctic snow samples, Nature, 184, 1557–1558, 1959. a

Gonfiantini, R.: Some results on oxygen isotope stratigraphy in the deep drilling at King Baudouin Station, Antarctica, J. Geophys. Res., 70, 1815–1819, 1965. a

Hartmann, D. L.: Global physical climatology, International Geophysics, Elsevier Science, ISBN 9780080918624, 2015. a

Hellmann, R. and Harvey, A. H.: First-principles diffusivity ratios for kinetic isotope fractionation of water in air, Geophys. Res. Lett., 47, e2020GL089999,, 2020. a, b, c

Hong, S.-Y., Dudhia, J., and Chen, S.-H.: A revised approach to ice microphysical processes for the bulk parameterization of clouds and precipitation, Mon. Weather Rev., 132, 103–120, 2004. a, b

Hu, Y., Rodier, S., Xu, K.-m., Sun, W., Huang, J., Lin, B., Zhai, P., and Josset, D.: Occurrence, liquid water content, and fraction of supercooled water clouds from combined CALIOP/IIR/MODIS measurements, J. Geophys. Res.-Atmos., 115, D00H34,, 2010. a, b, c, d, e

IAEA: GNIP Maps and Animations, International Atomic Energy Agency, Vienna, (last access: 13 May​​​​​​​ 2020), 2001. a, b, c

Johnsen, S. J., Clausen, H. B., Cuffey, K. M., Hoffmann, G., Schwander, J., and Creyts, T.: Diffusion of stable isotopes in polar firn and ice: the isotope effect in firn diffusion, Phys. Ice Core Rec., 159, 121–140, 2000. a

Jouzel, J. and Merlivat, L.: Deuterium and oxygen 18 in precipitation: Modeling of the isotopic effects during snow formation, J. Geophys. Res.-Atmos., 89, 11749–11757, 1984. a, b, c, d, e, f, g, h, i, j, k

Jouzel, J., Merlivat, L., and Lorius, C.: Deuterium excess in an East Antarctic ice core suggests higher relative humidity at the oceanic surface during the last glacial maximum, Nature, 299, 688–691, 1982. a, b

Jouzel, J., Alley, R. B., Cuffey, K., Dansgaard, W., Grootes, P., Hoffmann, G., Johnsen, S. J., Koster, R., Peel, D., Shuman, C., Stievenard, M., Stuiver, M., and White, J.: Validity of the temperature reconstruction from water isotopes in ice cores, J. Geophys. Res.-Oceans, 102, 26471–26487, 1997. a, b, c, d, e, f

Kahle, E. C., Steig, E. J., Jones, T. R., Fudge, T., Koutnik, M. R., Morris, V. A., Vaughn, B. H., Schauer, A. J., Stevens, C. M., Conway, H., Waddington, E. D., Buizert, C., Epifanio, J., and White, J. W. C.: Reconstruction of temperature, accumulation rate, and layer thinning from an ice core at South Pole, using a statistical inverse method, J. Geophys. Res.-Atmos., 126, e2020JD033300,, 2021. a, b, c, d

Kalnay, E., Kanamitsu, M., Kistler, R., Collins, W., Deaven, D., Gandin, L., Iredell, M., Saha, S., White, G., Woollen, J., Zhu, Y., Chelliah, M., Ebisuzaki, W., Higgins, W., Janowiak, J., Mo, K. C., Ropelewski, C., Wang, J., Leetmaa, A., Reynolds, R., Jenne, R., and Joseph, D.​​​​​​​: The NCEP/NCAR 40-year reanalysis project, B. Am. Meteorol. Soc., 77, 437–471, 1996. a, b

Kavanaugh, J. L. and Cuffey, K. M.: Generalized view of source region effects on δD and deuterium excess of ice sheet precipitation, Ann. Glaciol., 35, 111–117, 2002. a, b, c, d

Kavanaugh, J. L. and Cuffey, K. M.: Space and time variation of δ18O and δD in Antarctic precipitation revisited, Global Biogeochem. Cy., 17, 1017,, 2003. a, b, c, d, e, f, g, h

Kurita, N., Hirasawa, N., Koga, S., Matsushita, J., Steen-Larsen, H. C., Masson-Delmotte, V., and Fujiyoshi, Y.: Influence of large-scale atmospheric circulation on marine air intrusion toward the East Antarctic coast, Geophys. Res. Lett., 43, 9298–9305, 2016. a, b, c

Lamb, K. D., Clouser, B. W., Bolot, M., Sarkozy, L., Ebert, V., Saathoff, H., Möhler, O., and Moyer, E. J.: Laboratory measurements of HDO/H2O isotopic fractionation during ice deposition in simulated cirrus clouds, P. Natl. Acad. Sci. USA, 114, 5612–5617,, 2017. a

Landais, A., Barkan, E., and Luz, B.: Record of δ18O and 17O-excess in ice from Vostok Antarctica during the last 150,000 years, Geophys. Res. Lett., 35, L02709, 2008. a

Langway, C. C.: A 400 Meter Deep Ice Core in Greenland: Preliminary Report, J. Glaciol., 3, 217–217, 1958. a

Liu, J., Xiao, C., Ding, M., and Ren, J.: Variations in stable hydrogen and oxygen isotopes in atmospheric water vapor in the marine boundary layer across a wide latitude range, J. Environ. Sci., 26, 2266–2276, 2014. a, b, c, d, e, f

Luz, B., Barkan, E., Yam, R., and Shemesh, A.: Fractionation of oxygen and hydrogen isotopes in evaporating water, Geochim. Cosmochim. Acta, 73, 6697–6703, 2009. a, b, c, d, e, f, g

Majoube, M.: Fractionation factor of 18 O between water vapour and ice, Nature, 226, 1242–1242, 1970. a

Majoube, M.: Fractionnement en oxygène 18 et en deutèrium entre l'eau et sa vapeur, J. Chim. Phys., 68, 1423–1436, 1971. a

MARGO Project Members: Constraints on the magnitude and patterns of ocean cooling at the Last Glacial Maximum, Nat. Geosci., 2, 127–132,, 2009. a

Markle, B.: Simple Water Isotope Model, GitHub [code],, last access: 3 May​​​​​​​ 2022. a

Markle, B. R. and Steig, E. J.: bradley-markle/simple_water_isotope_model: Simple Water Isotope Model (Version v0), Zenodo [data set],, 2022. a

Markle, B., Bertler, N., Sinclair, K., and Sneed, S.: Synoptic variability in the Ross Sea region, Antarctica, as seen from back-trajectory modeling and ice core analysis, J. Geophys. Res.-Atmos., 117, D02113,, 2012. a

Markle, B. R., Steig, E. J., Buizert, C., Schoenemann, S. W., Bitz, C. M., Fudge, T., Pedro, J. B., Ding, Q., Jones, T. R., White, J. W., and Sowers, T.: Global atmospheric teleconnections during Dansgaard-Oeschger events, Nat. Geosci., 10, 36–40, 2017. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s, t, u, v, w, x

Markle, B. R., Steig, E. J., Roe, G. H., Winckler, G., and McConnell, J. R.: Concomitant variability in high-latitude aerosols, water isotopes and the hydrologic cycle, Nat. Geosci., 11, 853–859,, 2018. a, b

Masson-Delmotte, V., Kageyama, M., Braconnot, P., Charbit, S., Krinner, G., Ritz, C., Guilyardi, E., Jouzel, J., Abe-Ouchi, A., Crucifix, M., Gladstone, R. M., Hewitt, C. D., Kitoh, A., LeGrande, A. N., Marti, O., Merkel, U., Motoi, T., Ohgaito, R., Otto-Bliesner, B., Peltier, W. R., Ross, I., Valdes, P. J., Vettoretti, G., Weber, S. L., Wolk, F., and Yu, Y.: Past and future polar amplification of climate change: climate model intercomparisons and ice-core constraints, Clim. Dynam., 26, 513–529, 2006. a, b

Masson-Delmotte, V., Hou, S., Ekaykin, A., Jouzel, J., Aristarain, A., Bernardo, R., Bromwich, D., Cattani, O., Delmotte, M., Falourd, S., Frezzotti, M., Gallée, H., Genoni, L., Isaksson, E., Landais, A., Helsen, M. M., Hoffmann, G., Lopez, J., Morgan, V., Motoyama, H., Noone, D., Oerter, H., Petit, J. R., Royer, A., Uemura, R., Schmidt, G. A., Schlosser, E., Simões, J. C., Steig, E. J., Stenni, B., Stievenard, M., van den Broeke, M. R., van de Wal, R. S. W., van de Berg, W. J., Vimeux, F., and White, J. W. C.: A review of Antarctic surface snow isotopic composition: Observations, atmospheric circulation, and isotopic modeling, J. Climate, 21, 3359–3387, 2008. a, b, c, d, e, f, g, h, i, j

Merlivat, L.: Molecular diffusivities of H216O, HD16O, and H218O in gases, The J. Chem. Phys., 69, 2864–2871, 1978. a, b, c

Merlivat, L. and Jouzel, J.: Global climatic interpretation of the deuterium-oxygen 18 relationship for precipitation, J. Geophys. Res.-Oceans, 84, 5029–5033, 1979. a, b, c, d, e, f, g, h, i, j

Merlivat, L. and Nief, G.: Fractionnement isotopique lors des changements d ètat solide-vapeur et liquide-vapeur de l'eau à des tempèratures infèrieures à 0 C, Tellus, 19, 122–127, 1967. a

Murphy, D. and Koop, T.: Review of the vapour pressures of ice and supercooled water for atmospheric applications, Q. J. Roy. Meteorol. Soc., 131, 1539–1565, 2005. a

Murray, F. W.: On the computation of saturation vapor pressure, J. Appl. Meteorol. Clim., 6, 203–204,<0203:OTCOSV>2.0.CO;2, 1966. a

Otto-Bliesner, B. L., Brady, E. C., Clauzet, G., Tomas, R., Levis, S., and Kothavala, Z.: Last glacial maximum and Holocene climate in CCSM3, J. Climate, 19, 2526–2544, 2006. a

Pang, H., Hou, S., Landais, A., Masson-Delmotte, V., Prie, F., Steen-Larsen, H. C., Risi, C., Li, Y., Jouzel, J., Wang, Y., He, J., Minster, B., and Falourd, S.: Spatial distribution of 17O-excess in surface snow along a traverse from Zhongshan station to Dome A, East Antarctica, Earth Planet. Sci. Lett., 414, 126–133, 2015. a

Petit, J.-R., Jouzel, J., Raynaud, D., Barkov, N. I., Barnola, J.-M., Basile, I., Bender, M., Chappellaz, J., Davis, M., Delaygue, G., Delmotte, M., Kotlyakov, V. M., Legrand, M., Lipenkov, V. Y., Lorius, C., Pépin, L., Ritz, C., Saltzman, E., and Stievenard, M.: Climate and atmospheric history of the past 420,000 years from the Vostok ice core, Antarctica, Nature, 399, 429–436, 1999. a, b, c, d

Pfahl, S. and Wernli, H.: Lagrangian simulations of stable isotopes in water vapor: An evaluation of nonequilibrium fractionation in the Craig-Gordon model, J. Geophys. Res.-Atmos., 114, D20108,, 2009. a

Rayleigh, L.: LIX. On the distillation of binary mixtures, The London, Edinburgh, and Dublin Philosophical Magazine and Journal of Science, 4, 521–537, 1902. a

Risi, C., Landais, A., Bony, S., Jouzel, J., Masson-Delmotte, V., and Vimeux, F.: Understanding the 17O excess glacial-interglacial variations in Vostok precipitation, J. Geophys. Res.-Atmos., 115, D10112,, 2010. a, b, c, d, e, f

Schilla, A. S.: The stable isotopes and deuterium excess from the Siple Dome ice core: implications for the late Quaternary climate and elevation history of the Ross Sea Region, West Antarctica, Ph.D. thesis, University of Colorado at Boulder, ISBN 978-0-549-14257-7, 2007. a

Schmidt, G., Bigg, G., and Rohling, E.: Global seawater oxygen-18 database–v1. 21, http://data. giss. nasa. gov/o18data (last access: October 2016), 1999. a, b, c, d

Schoenemann, S. W., Steig, E. J., Ding, Q., Markle, B. R., and Schauer, A. J.: Triple water-isotopologue record from WAIS Divide, Antarctica: Controls on glacial-interglacial changes in 17Oexcess of precipitation, J. Geophys. Res.-Atmos., 119, 8741–8763, 2014. a, b, c, d

Sodemann, H. and Stohl, A.: Asymmetries in the moisture origin of Antarctic precipitation, Geophys. Res. Lett., 36, L22803,, 2009. a, b

Steig, E. J., Ding, Q., White, J. W., Küttel, M., Rupper, S. B., Neumann, T. A., Neff, P. D., Gallant, A. J., Mayewski, P. A., Taylor, K. C., Hoffmann, G., Dixon, D. A., Schoenemann, S. W., Markle, B. R., Fudge, T. J., Schneider, D. P., Schauer, A. J., Teel, R. P., Vaughn, B. H., Burgener, L., Williams, J., and Korotkikh, E.: Recent climate and ice-sheet changes in West Antarctica compared with the past 2,000 years, Nat. Geosci., 6, 372–375, 2013. a

Steig, E. J., Jones, T. R., Schauer, A. J., Kahle, E. C., Morris, V. R., Vaughn, B. H., Davidge, L., and White, J. W. C.: Continuous-flow analysis of δ17O, δ18O, and δD of H2O on an ice core from the South Pole, Front. Earth Sci., 9, 640292,, 2021. a

Stenni, B., Jouzel, J., Masson-Delmotte, V., Röthlisberger, R., Castellano, E., Cattani, O., Falourd, S., Johnsen, S., Longinelli, A., Sachs, J., Selmo, E., Souchez, R., Steffensen, J. P., and Udist, R.i: A late-glacial high-resolution site and source temperature record derived from the EPICA Dome C isotope records (East Antarctica), Earth Planet. Sci. Lett., 217, 183–195, 2004. a, b, c, d

Stenni, B., Masson-Delmotte, V., Selmo, E., Oerter, H., Meyer, H., Röthlisberger, R., Jouzel, J., Cattani, O., Falourd, S., Fischer, H., Hoffmann, G., Iacumin, P., Johnsen, S. J., Minster, B., and Udisti, R.: The deuterium excess records of EPICA Dome C and Dronning Maud Land ice cores (East Antarctica), Quaternary Sci. Rev., 29, 146–159, 2010. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o

Stenni, B., Buiron, D., Frezzotti, M., Albani, S., Barbante, C., Bard, E., Barnola, J., Baroni, M., Baumgartner, M., Bonazza, M., Capron, E., Castellano, E., Chappellaz, J., Delmonte, B., Falourd, S., Genoni, L., Iacumin, P., Jouzel, J., Kipfstuhl, S., Landais, A., Lemieux-Dudon, B., Maggi, V., Masson-Delmotte, V., Mazzola, C., Minster, B., Montagnat, M., Mulvaney, R., Narcisi, B., Oerter, H., Parrenin, F., Petit, J. R., Ritz, C., Scarchilli, C., Schilt, A., Schüpbach, S., Schwander, J., Selmo, E., Severi, M., Stocker, T. F., and Udisti, R.: Expression of the bipolar see-saw in Antarctic climate records during the last deglaciation, Nat. Geosci., 4, 46–49, 2011. a, b, c

Tegen, I. and Fung, I.: Modeling of mineral dust in the atmosphere: Sources, transport, and optical thickness, J. Geophys. Res.-Atmos., 99, 22897–22914, 1994.  a

Uemura, R., Matsui, Y., Yoshimura, K., Motoyama, H., and Yoshida, N.: Evidence of deuterium excess in water vapor as an indicator of ocean surface conditions, J. Geophys. Res.-Atmos., 113, D19114,, 2008. a, b, c, d, e, f, g, h

Uemura, R., Barkan, E., Abe, O., and Luz, B.: Triple isotope composition of oxygen in atmospheric water vapor, Geophys. Res. Lett., 37, L04402,, 2010. a, b, c, d, e, f, g

Uemura, R., Masson-Delmotte, V., Jouzel, J., Landais, A., Motoyama, H., and Stenni, B.: Ranges of moisture-source temperature estimated from Antarctic ice cores stable isotope records over glacial–interglacial cycles, Clim. Past, 8, 1109–1125,, 2012. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q, r, s, t, u, v, w, x, y, z, aa, ab, ac, ad, ae

Vimeux, F., Cuffey, K. M., and Jouzel, J.: New insights into Southern Hemisphere temperature changes from Vostok ice cores using deuterium excess correction, Earth Planet. Sci. Lett., 203, 829–843,, 2002. a, b, c, d, e, f, g, h

WAIS Divide Project Members: Onset of deglacial warming in West Antarctica driven by local orbital forcing, Nature, 500, 440–444,, 2013. a, b, c

WAIS Divide Project Members: Precise interpolar phasing of abrupt climate change during the last ice age, Nature, 520, 661–665,, 2015. a

Werner, M., Jouzel, J., Masson-Delmotte, V., and Lohmann, G.: Reconciling glacial Antarctic water stable isotopes with ice sheet topography and the isotopic paleothermometer, Nat. Commun., 9, 1–10, 2018. a, b

Winkler, R., Landais, A., Sodemann, H., Dümbgen, L., Prié, F., Masson-Delmotte, V., Stenni, B., and Jouzel, J.: Deglaciation records of 17O-excess in East Antarctica: reliable reconstruction of oceanic normalized relative humidity from coastal sites, Clim. Past, 8, 1–16,, 2012. a

Xiao, C., Ding, M., Masson-Delmotte, V., Zhang, R., Jin, B., Ren, J., Li, C., Werner, M., Wang, Y., Cui, X., and Wang, X.: Stable isotopes in surface snow along a traverse route from Zhongshan station to Dome A, East Antarctica, Clim. Dynam., 41, 2427–2438, 2013. a


This process is not incidental but rather fundamental to the climate system itself. The Earth's surface absorbs shortwave radiation from the sun and transfers that heat to the atmosphere. The majority of that transferred heat is latent in the form of evaporated moisture. The basic climatological function of the atmosphere and its motions is to transport heat from regions of net energy gain to regions of net loss from low to high latitudes, and a large fraction of that transported heat is also moisture.


The base model includes local closure during evaporation, values of SST0 and RH0 for a given T0 determined by a spline fit to NCEP/NCAR climatology, and a tuned supersaturation such that S=1-0.00525T, where T is in degrees Celsius (C). See the Appendix for details.

Short summary
The geochemistry preserved in polar ice can provide detailed histories of Earth’s climate over millennia. Here we use the stable isotope ratios of ice from many Antarctic ice cores to reconstruct temperature variability of Antarctica and the midlatitude Southern Hemisphere over tens of thousands of years. We improve upon existing methods to estimate temperature from the geochemical measurements and investigate the patterns of climate change in the past.