the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
Improving temperature reconstructions from icecore waterisotope records
Bradley R. Markle
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 icecore records generally rely on linear approximations of the relationships among local temperature, source temperature, and waterisotope values. However, there are important nonlinearities that significantly affect such reconstructions, particularly for source region temperatures. Here, we describe a relatively simple waterisotope 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 waterisotope 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 icecore 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.
Stableisotope ratios of water have been the foundational proxy for polar paleoclimate research for more than a halfcentury (Langway, 1958; Gonfiantini, 1959; Dansgaard, 1964). Primarily used as a temperature proxy, stratigraphic records of waterisotope 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 (MassonDelmotte et al., 2006; EPICA Community Members, 2006; WAIS Divide Project Members, 2013, 2015). Both oxygen and hydrogen have stable isotopes whose ratios (${}^{\mathrm{18}}\mathrm{O}/{}^{\mathrm{16}}\mathrm{O}$ and ${}^{\mathrm{2}}\mathrm{H}/{}^{\mathrm{1}}\mathrm{H}$) are commonly expressed as deviations, δ^{18}O and δD, from Vienna Standard Mean Ocean Water (VSMOW) in per mill (‰):
where R_{x} is the ratio in the sample and R_{std} 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 waterisotope ratios during phase changes are all processes inherently linked to temperature and together underpin the use of waterisotope ratios in polar precipitation as a temperature proxy (Craig, 1961; Epstein et al., 1963; Dansgaard, 1964; Gonfiantini, 1965). The strong empirical correlation between the waterisotope ratios in precipitation and surface temperature supports this interpretation (Petit et al., 1999; Jouzel et al., 1997; MassonDelmotte 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 icecore waterisotope records if the relevant scaling relationships can be determined from theory, models, or observations (Vimeux et al., 2002; Kavanaugh and Cuffey, 2002; Stenni et al., 2010). Here, we examine the widely used assumption of linearity in the scaling relationships between waterisotope ratios and temperature.
1.1 Temperature reconstructions
Any interpretation of waterisotope 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., Dansgaard, 1964) is sufficient to qualitatively interpret variations in waterisotope ratios as variations in temperature in the high latitudes. The simplest quantitative interpretation of icecore waterisotope records relies on the empirical correlation between observed waterisotope 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 highlatitude waterisotope 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 subequilibrium relative humidity and influenced by sea surface temperature and wind speed (Merlivat and Jouzel, 1979; Jouzel et al., 1982).
A more complete approach to reconstructing temperature from waterisotope 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 δ^{18}O and δD, d_{xs}=δD$\mathrm{8}\times {\mathit{\delta}}^{\mathrm{18}}$O, and is commonly used to quantify these differences (Dansgaard, 1964; Merlivat and Jouzel, 1979).
Changes in waterisotope parameters measured in precipitation at an icecore site, Δδ^{18}O and Δd_{xs}, can be conceptualized as driven by changes in site and evaporation source temperature, ΔT_{site} and ΔT_{source}:
where β and γ are the partial derivatives of δ^{18}O and d_{xs} with site and source temperature, respectively. The magnitudes of β and γ can be diagnosed from waterisotope distillation models for the icecore site in question (Vimeux et al., 2002; Kavanaugh and Cuffey, 2002; Stenni et al., 2010; Uemura et al., 2012). Once these slopes are established, the equations may be solved for ΔT_{site} and ΔT_{source} using records of δ^{18}O and d_{xs} (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 icecore 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 (Dansgaard, 1964). However, a linear relationship between δ^{18}O and δD is not fundamental (Craig, 1961); equilibrium fractionation alone drives a nonlinear relationship between δ^{18}O 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 Jouzel, 1979; 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 d_{xs} 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 T_{site}.
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 secondorder polynomial to a compilation of $\mathit{\delta}{{}^{\prime}}^{\mathrm{18}}\mathrm{O}$ and δ^{′}D data, where ${\mathit{\delta}}_{x}{}^{\prime}=\mathrm{ln}\left(\mathrm{1}+{\mathit{\delta}}_{x}\right)$, and defined a phenomenological, nonlinear deuterium excess parameter:
with coefficients $A=\mathrm{28.5}$ and B=8.47 (note that the coefficients and δ^{′} values are unitless; for example, $\mathit{\delta}{{}^{\prime}}^{\mathrm{18}}\mathrm{O}=\mathrm{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 d_{xs} 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 waterisotopebased temperature reconstructions and suggest an improved technique.
The quantitative reconstruction of temperatures from waterisotope 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 waterisotope 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 highlatitude waterisotope 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 longterm mean is strongly influenced by climatological moisture distillation.^{1}
Our model distills moisture down thermodynamic pathways defined by temperature and pressure. Changes in waterisotope 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, T_{0}, to a final condensation temperature, T_{c}. The pathways are pseudoadiabatic, 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, T_{0}, to other initial conditions including sea surface temperature, SST_{0}, and relative humidity, RH_{0} (see Sect. A1.1).
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 pseudoadiabatic pathway, initial climatological correlations during evaporation, and nonuniqueness. 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, precipitationweighted condensation temperature above that site. We make no attempt to model postdepositional processes.
To investigate the sensitivity of waterisotope 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={T}_{\mathrm{0}},{T}_{\mathrm{0}}\mathrm{d}T,\mathrm{\dots},{T}_{\mathrm{c}}$, with dT=0.1 ^{∘}C. We run nearly 24 000 trajectories to fill the plausible parameter space of 0 ^{∘}C $\le {T}_{\mathrm{0}}\le \mathrm{28}$ ^{∘}C and −70 ^{∘}C $\le {T}_{\mathrm{c}}\le \mathrm{10}$ ^{∘}C. We first examine the model with base assumptions and parameterizations^{2}. 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, T_{c} and T_{0}, respectively, and whose z (color) dimension is the isotopic value of final precipitation in Antarctica: δ^{18}O in Fig. 3a, δD in Fig. 3b, d_{xs} in Fig. 3c, and d_{ln} in Fig. 3d.
The gradients of both the δ^{18}O 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 δ^{18}O and δD are not strictly linear with condensation temperature T_{c}, clearly varying with its absolute value (and to a much lesser extent with the evaporation temperature, T_{0}, due to its influence on the total distillation gradient). Further, the partial slopes of δ^{18}O 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 T_{0} and T_{c} are shown explicitly in Figs. 4 and 5. It is important to recognize that the partial derivatives with respect to T_{0} 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.
The modeled d_{xs} surface shows strong slopes along both the condensation temperature and evaporation temperature axes (Fig. 3c), as does modeled d_{ln} (Fig. 3d). The d_{ln} depends more strongly on the evaporation temperature than the d_{xs}. In particular, at the coldest condensation temperatures, variability in d_{xs} 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 (d_{ln}) is a more faithful qualitative proxy for source region conditions than the linear definition (d_{xs}). Even at very low condensation temperatures, d_{ln} still depends strongly on the initial evaporation temperature, whereas linear d_{xs} 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 −30 ^{∘}C, 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 T_{c} is smoothed somewhat when atmospheric mixing is incorporated into the model (see Appendix A5). The changes in the partial derivatives with respect to T_{c} across the entire range investigated are larger than these changes localized around −30 ^{∘}C.
The temperature reconstruction technique described in Sect. 1.1 is based on the linearization of the slopes between δ^{18}O, d_{xs}, condensation temperature, and evaporation temperature, and it assumes that the β and γ parameters in Eqs. (2) and (3) are fixed over the range of reconstructed ΔT_{site} and ΔT_{source}. 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 $\frac{\partial {\mathit{\delta}}^{\mathrm{18}}\mathrm{O}}{\partial {T}_{\mathrm{c}}}$ and $\frac{\partial {\mathit{\delta}}^{\mathrm{18}}\mathrm{O}}{\partial {T}_{\mathrm{0}}}$ in Figs. 4 and 5. Although the slope of δ^{18}O along the condensation temperature axis, $\frac{\partial {\mathit{\delta}}^{\mathrm{18}}\mathrm{O}}{\partial {T}_{\mathrm{c}}}$, does not change dramatically, it is clearly variable, as is $\frac{\partial {\mathit{\delta}}^{\mathrm{18}}\mathrm{O}}{\partial {T}_{\mathrm{0}}}$. The slopes $\frac{\partial {d}_{\mathrm{xs}}}{\partial {T}_{\mathrm{c}}}$ and $\frac{\partial {d}_{\mathrm{xs}}}{\partial {T}_{\mathrm{0}}}$ (comparable to β_{1} and β_{2} in Eq. 3, respectively) are highly variable. Indeed, the d_{xs} surface in Fig. 3 has a saddle at moderate condensation temperatures, over which $\frac{\partial {d}_{\mathrm{xs}}}{\partial {T}_{\mathrm{c}}}$ changes sign (Fig. 4c). This shows that the assumption of constant β and γ parameters in isotopebased temperature reconstructions is incorrect and may be reasonable only under narrow ranges of ΔT_{site} and ΔT_{source}. 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 T_{0}, the evaporation source temperature.
Evidence for the critical change in the sign of $\frac{\partial {d}_{\mathrm{xs}}}{\partial {T}_{\mathrm{c}}}$ (Fig. 4c) can actually be seen directly in Antarctic icecore records. Deuterium excess (d_{xs}) 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 d_{xs} 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, d_{ln} values from different cores are vastly more coherent (Fig. 6e and f) for reasons described above. Comparison of d_{xs} and d_{ln} between cores (Figs. 6 and A24) together with our modeling results demonstrates that the incoherence of Antarctic d_{xs} records arises from the change in slope of $\frac{\partial {d}_{\mathrm{xs}}}{\partial {T}_{\mathrm{c}}}$, not from differing source conditions between sites. This exposes a fundamental flaw in the assumptions used in traditional isotope temperature reconstructions.
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 icecore 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 waterisotope–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 waterisotope fractionation process.
4.1 Reconstruction method
For every pair of T_{0} and T_{c} inputs to SWIM there is a corresponding modeled value of δ^{18}O, δD, and d_{ln} 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., δ^{18}O and d_{ln}. This defines a set of maps, with x and y axes of δ^{18}O and d_{ln} and with z axes T_{c} and T_{0}, as shown in Fig. 7. To reconstruct T_{c} and T_{0}, the inverted model results may be used as a lookup table: a pair of δ^{18}O and d_{ln} measurements determine a pair of T_{c} and T_{0} reconstructions (Fig. A16). While previous reconstruction methods (e.g., Vimeux et al., 2002; Kavanaugh and Cuffey, 2002; Stenni et al., 2010) linearize the slopes calculated by a waterisotope 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 highorder polynomials to the results.
The boundary conditions T_{c} and T_{0} 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 δ^{18}O, δD measurements (d_{xs} and d_{ln} being secondorder parameters), any combination of those parameters may seem equally well suited for the purposes of temperature reconstruction. In practice, however, δ^{18}O and d_{ln} 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 d_{ln} is a better qualitative proxy for source region temperature than d_{xs}. After proposing the d_{ln} parameter, Uemura et al. (2012) suggest that there is no added value in the logarithmic parameter over the traditional linear d_{xs} 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 (RH_{0}) instead of source region air temperature (T_{0}), since we assume climatological relationships between them. Indeed, RH_{0} is a strong lever on the deuterium excess of evaporation. We reconstruct T_{0} 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 T_{c}, T_{0}, and RH_{0}, by modeling multidimensional parameter spaces. This of course requires additional measured constraints, which could include the δ^{17}O 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 highresolution MERRA2 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 temperaturedependent 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 nonuniqueness on our temperature reconstruction technique arising from belowfreezing 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.
We reconstruct condensation site and surface temperatures and evaporation source temperatures for eight different Antarctic deep icecore sites for which there are δ^{18}O and d_{ln} records (Fig. 9). The records include WDC (Markle et al., 2017; WAIS Divide Project Members, 2013; Steig et al., 2013) and Siple Dome (Brook et al., 2005; Schilla, 2007) 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 Wal, 2008), δ^{18}O_{sw}, 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 Members, 2013). 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 waterisotope 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 icecore sites by regressing the T_{c} and T_{0} temperature fields to subsets of δ^{18}O and d_{xs} 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 icecore 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 icecore 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 ±1 ^{∘}C 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 ±2 ^{∘}C. Further, the total variability in reconstructed evaporation temperature is much smaller than that in icecore 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 T_{s} and especially of T_{0} (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 d_{ln} parameter over the linear d_{xs} parameter (Fig. 6; Markle et al., 2017).
The relationship between δ^{18}O and T_{c} is largely linear across a wide range of values of T_{c}, regardless of evaporation temperature. The icecore 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 δ^{18}OtoT_{c} 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 liquidonly and iceonly condensate. SWIM retains liquid condensate at colder temperatures than previous models (e.g., Kavanaugh and Cuffey, 2002), in line with satellite measurements (Hu et al., 2010). The resulting transition of fractionation factors drives the nonlinearities in the δ^{18}OtoT_{c} 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 δ^{18}O and T_{c} 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 waterisotope 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.
Using the selfconsistent, nonlinear temperature reconstruction technique for eight different icecore 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 icecore 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.8 ^{∘}C warming at WDC during the deglaciation; we reconstruct 11.2±0.5 ^{∘}C 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.68 ^{∘}C between 19.5–22.5 and 0.5–2.5 ka. Our reconstruction yields a site temperature warming of 8.9±0.4 ^{∘}C 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.
We create a stack of each reconstructed temperature variable (the evaporation temperature T_{0}, condensation temperature T_{c}, and surface temperature T_{s}) for all eight icecore 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 Antarcticwide 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.
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 icecore records as a function of the mean latitude of the moisture source distribution for each site (based on a watertagged 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 T_{0}–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 icecore site. The moisture source points are plotted at the longitude of the respective icecore 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 T_{0} 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 icecore T_{0} reconstructions have low spatial resolution owing to broad moisture source distributions, they benefit from the temporal resolution and precision of the icecore 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 (MassonDelmotte et al., 2006; OttoBliesner 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 Members, 2009, 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 multimodel 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 icecore 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., >20 ^{∘}C), 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 GCMmodeled changes in Antarctic water isotopes and icecore records. Understanding the full set of processes responsible for the reconstructed pattern of Southern Hemisphere polar amplification is a topic for future work.
Icecore 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, δ^{18}O 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, d_{xs}, is an unreliable indicator of relative evaporation site temperature change, particularly at East Antarctic sites with very depleted δ^{18}O and δD values. In these cases, the logarithmic definition of the parameter, d_{ln}, is a more faithful qualitative proxy for evaporation temperature.
We can use models to make quantitative interpretations of waterisotope variability and to disentangle the combined influences of the source and site temperatures. To date, most waterisotope temperature inversions have assumed linear relationships (Kavanaugh and Cuffey, 2003; Vimeux et al., 2002; Stenni et al., 2011; Uemura et al., 2012). However, as shown here, this assumption is flawed. Even in the simplified waterisotope 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 longstanding debate regarding the interpretation of “spatial” and “temporal” slopes in the waterisotope–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 waterisotope fractionation. Neither space nor time can independently cause water to change phase and fractionate.
The fundamental dimension through which to understand waterisotope 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 waterisotope–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 waterisotope 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.
The Simple Water Isotope Model (SWIM) is based on existing numerical Rayleightype distillation models (Merlivat and Jouzel, 1979; Jouzel and Merlivat, 1984; Ciais and Jouzel, 1994; Criss, 1999; Kavanaugh and Cuffey, 2003), 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 pseudoadiabatic temperature pathways from an initial surface air temperature, T_{0}, to a final condensation temperature, T_{c}, and discrete steps dT, as well as Euler numerics.
A1.1 Source region conditions
The moisture source surface air temperature (T_{a}), 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 SST_{0} and RH_{0} given a specified initial air temperature, T_{0}, using the 1980–2010 annual mean climatological fields from the NCEP/NCAR reanalysis project (Kalnay et al., 1996) and the ERAInterim 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 welldefined 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 (Dai, 2006; Vimeux et al., 2002). We find that the overocean 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, T_{0}, our model uses values of SST_{0} and RH_{0} based on fits to the modern climatology. We use three methods to calculate the climatological relationships over the interval −10^{∘} $\le {T}_{\mathrm{a}}\le $ 28 ^{∘}C: a cubic spline with specified noise tolerance, the mean and median of SST and RH distributions within binned values of T_{a}, and highorder 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 T_{a}toSST fit between the Northern and Southern Hemisphere for air temperatures between 5 and −15 ^{∘}C, as seen in Fig. A1. In this study we will use the Southern Hemisphere fit for T_{a} to SST. We find no major hemispheric differences for the T_{a}toRH fit and find little impact of zonal asymmetry on either fit. We find relatively small differences in the fit between T_{a} and SST_{0} for different seasons and somewhat larger seasonal changes in the T_{a} and RH_{0} relationship. We use the annual average fits hereafter and test the sensitivity of the waterisotope values of evaporation to these seasonal differences.
The normalized relative humidity, RH_{n}, is critical to kinetic fractionation of water isotopes during evaporation (Merlivat and Jouzel, 1979; Risi et al., 2010) and depends on the three variables above:
where e_{s} (T_{a}) and e_{s} (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, T_{0}, and specified surface pressure, P_{0}, moisture is transported toward the pole in isolation, cooling and condensing along the way. The air parcel is cooled pseudoadiabatically, 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 pseudoadiabatic 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 pseudoadiabat 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 temperaturedependent saturated vapor pressures of ice and liquid water (Murray, 1966; Murphy and Koop, 2005), together with air pressure P(T), define saturated mixing ratios for ice and liquid water,
as functions of T, where $\frac{{R}_{\mathrm{d}}}{{R}_{\mathrm{wv}}}$ is the ratio of gas constants of dry air and water vapor.
We consider air parcels with mixed ice and liquid condensate (Ciais and Jouzel, 1994), in which the ice fraction smoothly increases as temperatures decreases below freezing. Many models, including isotopeenabled 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 Jouzel, 1994; Kavanaugh and Cuffey, 2003). We use temperaturedependent 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 Cuffey, 2003).
The specific heat at constant pressure, c_{p}, latent heat, L, saturated vapor pressure, e_{s}, and saturated mixing ratio, r_{s}, are temperaturedependent 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 Cuffey, 2003). For example, ${r}_{\mathrm{s}\left(\mathrm{eff}\right)}={r}_{\mathrm{s}\left(\mathrm{ice}\right)}{F}_{\left(\mathrm{ice}\right)}+{r}_{\mathrm{s}\left(\mathrm{liq}\right)}{F}_{\left(\mathrm{liq}\right)}$, where r_{s(ice)} and r_{s(liq)} are the saturated mixing ratios of ice and liquid, respectively, and F_{(ice)} and F_{(liq)} are the (temperaturedependent) fractions of each phase of condensate.
Moisture is removed along the temperature pathway owing to the temperaturedependent changes in the saturated mixing ratio, $\frac{\mathrm{d}q}{\mathrm{d}T}=\frac{\mathrm{d}{r}_{\mathrm{s}\left(\mathrm{eff}\right)}}{\mathrm{d}T}$. This is a simplified view of largescale 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 Merlivat, 1984). A paucity of condensation nuclei may lead to further supersaturation at very cold temperatures (Tegen and Fung, 1994). 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 waterisotope fractionation scheme used in SWIM. We model equilibrium and the kinetic fractionation of the $\frac{{}^{\mathrm{2}}\mathrm{H}}{{}^{\mathrm{1}}\mathrm{H}}$, $\frac{{}^{\mathrm{18}}\mathrm{O}}{{}^{\mathrm{16}}\mathrm{O}}$, and $\frac{{}^{\mathrm{17}}\mathrm{O}}{{}^{\mathrm{16}}\mathrm{O}}$ ratios in water. We use conventional notation in which R is the number ratio of heavy to light isotopes of a species, for example ${}^{D}R=\frac{{}^{\mathrm{2}}\mathrm{H}}{{}^{\mathrm{1}}\mathrm{H}}$ and ${}^{\mathrm{18}}R=\frac{{}^{\mathrm{18}}\mathrm{O}}{{}^{\mathrm{16}}\mathrm{O}}$.
The fractionation factor is the ratio of R values between phases. For example, the fractionation factor between liquid and vapor phases for δ^{18}O is
We use the empirically determined, temperaturedependent 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)} (Majoube, 1970, 1971; Merlivat and Nief, 1967; Criss, 1999), 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 δ^{18}O of seawater (δ^{18}O_{sw}) 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 δ^{18}O_{sw} and δD_{sw} from a compilation of global seawater measurements (Schmidt et al., 1999) to find an initial δD_{sw} given the specified initial δ^{18}O_{sw} (which changes with mean climate) and investigate the sensitivity of the model to these initial conditions. We assume a ^{17}O_{xs} of seawater equal to 0, where ${}^{\mathrm{17}}{\mathrm{O}}_{\mathrm{xs}}=\mathit{\delta}{{}^{\prime}}^{\mathrm{17}}\mathrm{O}\mathrm{0.528}\times \mathit{\delta}{{}^{\prime}}^{\mathrm{18}}$O, 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 % (Hartmann, 2015). Because of this steadystate 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 Jouzel, 1979). The effective fractionation factor associated with diffusion and turbulence is
where D and D^{*} are the diffusivities of the light and heavy isotopes, respectively (Merlivat and Jouzel, 1979). 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 ${\mathrm{H}}_{\mathrm{2}}^{\mathrm{18}}\mathrm{O}$ and ${\mathrm{H}}_{\mathrm{2}}^{\mathrm{16}}\mathrm{O}$ during initial evaporation, the fractionation factor ^{18}α_{diff} equals 1.0 for pure turbulence and 1.0028 for pure molecular diffusion (Merlivat and Jouzel, 1979; Barkan and Luz, 2007).
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, δ^{18}O, and δ^{17}O in vapor above the Southern Ocean. Uemura et al. (2010) find a value of ${}^{\mathrm{18}}{\mathit{\alpha}}_{\mathrm{diff}}=\mathrm{1.007}\pm \mathrm{0.0013}$ and 1.008±0.0018 when optimizing for observations of d_{xs} and ^{17}O_{xs} 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 d_{xs} and ^{17}O_{xs} of vapor when using the observed values of T_{0}, 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 (Merlivat, 1978; Luz et al., 2009).
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) (<10 ^{∘}C), we use the measured value at 10 ^{∘}C (ϕ_{diff}=1.06). The differences in model results for evaporation with a temperaturedependent ϕ_{diff} and a constant ϕ_{diff}=0.88 are small (<1 ‰ for initial δD of vapor).
For the fractionation of $\frac{{}^{\mathrm{17}}\mathrm{O}}{{}^{\mathrm{16}}\mathrm{O}}$ we use the following relationships, which are backed by both theory and empirical observation (Barkan and Luz, 2005, 2007): ${}^{\mathrm{17}}{\mathit{\alpha}}_{\mathrm{eq}}{=}^{\mathrm{18}}{\mathit{\alpha}}_{\mathrm{eq}}^{\mathrm{0.529}}$ for vapor and liquid in equilibrium and ${}^{\mathrm{17}}{\mathit{\alpha}}_{\mathrm{diff}}{=}^{\mathrm{18}}{\mathit{\alpha}}_{\mathrm{diff}}^{\mathrm{0.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 Harvey, 2020), 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, RH_{n}, and the equilibrium and diffusive fractionation factors, α_{eq(l−v)} and α_{diff}. Following Criss (1999) and Luz et al. (2009),
where R_{o} and R_{v} are the isotopic ratios of the ocean water and the water vapor in the atmospheric boundary layer, respectively. R_{e} is the ratio of the evaporate, the net vapor lost to the atmosphere, which is a quantity that is not directly measurable (Criss, 1999).
If we assume that the only source of vapor to the boundary layer is the local evaporate, we may equate R_{v} and R_{e} and solve Eq. (A6) for R_{v} (Merlivat and Jouzel, 1979; Criss, 1999; Risi et al., 2010):
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 δ^{18}O and too low in d_{xs}, and these offsets are a function of environmental conditions (Risi et al., 2010).
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 δ^{18}O = −4.5 ‰ and δD $=\mathrm{26.7}$ ‰ (Craig and Gordon, 1965). Considering the global average steady state in which the ocean is the ultimate vapor source for precipitation (Merlivat and Jouzel, 1979; Criss, 1999) , the delta values of precipitation must reflect the net loss by evaporate from the ocean. Thus, globally, the ratio $\frac{{R}_{\mathrm{o}}}{{R}_{\mathrm{e}}}$ is 1.0267 for D and 1.0045 for ^{18}O. Instead of equating R_{e} and R_{v} locally, we define a global ${\stackrel{\mathrm{\u203e}}{\mathit{\alpha}}}_{\mathrm{evap}}={\left(\frac{{R}_{\mathrm{o}}}{{R}_{\mathrm{e}}}\right)}_{\mathrm{global}}$ and solve for R_{v}. Substitution into Eq. (A6) and rearranging leads to
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, R_{v}, is comprised of both local evaporate, R_{e}, and vapor evaporated at some distal location and advected to the site, ${R}_{\mathrm{v}}=(\mathrm{1}\mathit{\theta}){R}_{\mathrm{e}}+\mathit{\theta}{R}_{\mathrm{distal}}$, where θ is the fraction of nonlocal vapor with composition R_{distal}. The local evaporative conditions are defined by T_{0}, while the distal conditions may be either warmer or colder than T_{0}. Vapor is evaporated under those distal conditions (using a local closure assumption) and advected without fractionation before mixing with the local vapor of composition R_{e}. Figure A3 shows the isotopic composition of vapor over a range of T_{0}, with a full range of contributions from both a 5 ^{∘}C warmer moisture source (T_{0}+5 ^{∘}C, red) and a 5 ^{∘}C colder moisture source (T_{0}−5 ^{∘}C, 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.
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 wellmixed 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.
It is important to note that in Fig. A2 we use climatological correlations between T_{0}, SST, and RH, while the observational data represent far shorter time intervals, mostly from one season. When using the observed values of T_{0}, 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 d_{xs} of vapor and modeled ^{17}O_{xs} of vapor to Southern Ocean vapor observations using the observed environmental conditions at the time of the vapor measurements. The modeled relationships between d_{xs} and ^{17}O_{xs} with SST_{0} and RH_{0} 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 T_{0}, SST_{0}, and RH_{0}, 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 δ^{18}O_{sw} of the ocean (Fig. A5g–i) as well as the value of α_{diff} (Fig. A5j–l).
The direct comparisons of observed and modeled vapor composition using observed T_{0}, SST_{0}, and RH_{0} 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 temperaturedependent α_{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 shipmeasured RH and T_{a} and that felt at the water's surface or specific mixing of nonlocal moisture in the boundary layer on the days of the shipbased measurements. Given the magnitude of the misfit, it is also possible that spatial or temporal variability in δ^{18}O_{sw} and δD_{sw} 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 (Rayleigh, 1902; Dansgaard, 1964; Criss, 1999) is
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 r_{s} owing to pseudoadiabatic cooling from the source (Dansgaard, 1964) . Thus,
In general, condensation occurs in the model at saturation, and thus the temperaturedependent 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 Merlivat, 1984), the total fractionation α_{tot} factor is α_{tot}=α_{eq}α_{k}. Equation (A9) thus becomes
The kinetic fractionation factor, α_{k}, depends on the supersaturation of vapor over ice, S_{i}:
Following Jouzel and Merlivat (1984), we use the ratio of diffusivities for oxygen isotopes $\frac{{D}^{\mathrm{16}}}{{D}^{\mathrm{18}}}=\mathrm{1.0285}$ during moisture transport, representative of pure molecular diffusion and ignoring the negligible ventilation effect. Likewise we use $\frac{{D}^{\mathrm{1}}}{{D}^{\mathrm{2}}}=\mathrm{1.0251}$ for the ratio of diffusivities of hydrogen isotopes. These values imply a constant ϕ_{diff} during transport equal to 0.88 (Jouzel and Merlivat, 1984) rather than the temperaturedependent ϕ_{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 temperaturedependent 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 mixedphase 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),
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 waterisotope 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 waterisotope fractionation and the uncertainty in the underlying physics, the supersaturation is often parameterized to depend on temperature and tuned to fit waterisotope models to observations (Jouzel and Merlivat, 1984; Kavanaugh and Cuffey, 2003; Schoenemann et al., 2014). Jouzel and Merlivat (1984) parameterized the supersaturation as a function of temperature and note that available waterisotope data could not distinguish among possible functional forms of the parameterization (e.g., linear, exponential). Their linear parameterization has been used extensively in waterisotope models since:
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, iceonly conditions (Hong et al., 2004), rather than returning to unity as the cloud becomes entirely icephase. It is common for waterisotope 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 selfconsistent. Because the environmental supersaturation experienced by the air parcel is related to the relationship between temperature and moisture removal (that is $\frac{\mathrm{d}\mathrm{ln}\left(f\right)}{\mathrm{d}T}$) 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 waterisotope–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 r_{s(eff)}=r_{s(ice)}S_{i}. 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–waterisotope 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 Merlivat, 1984), 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 δ^{18}O and δD irreconcilable with observations. Were the air parcel to return to nonsupersaturated conditions in the iceonly portion of the cooling pathway, the simultaneous transition to equilibriumonly fractionation would drive a slope of $\frac{\partial {\mathit{\delta}}^{\mathrm{18}}\mathrm{O}}{\partial \mathit{\delta}\mathrm{D}}$ that is not compatible with measured values in Antarctic precipitation. This gives additional credence to sustained supersaturation at cold temperatures.
A3 Application to Antarctic icecore 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 moisturetagged 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 higherelevation sites having more equatorward moisture sources (Fig. 1d) due to topographic isolation (Sodemann and Stohl, 2009; Bailey et al., 2019). There are some notable asymmetries in this general pattern. The large embayments are areas of comparatively highlatitude 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 moisturetagged 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 Stohl, 2009; Markle et al., 2012; Buizert et al., 2018). These MSDs dictate the influence of evaporation source conditions (T_{a}, 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 icecore 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 icecore site, what is the meaning of the T_{source} that can be reconstructed from waterisotope records? It is the moistureweighted 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 waterisotope 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 waterisotope 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 T_{c} reconstructed from icecore 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 (Connolley, 1996). MassonDelmotte et al. (2008) review the relationship between the condensation temperature and the surface temperature (T_{s}) over Antarctica and compare the surface temperature to the weighted annual mean condensation temperature in both ERA40 reanalysis (1980–2002) and MAR, a highresolution mesoscale model forced by ERA40 (see Fig. 8, MassonDelmotte 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). MassonDelmotte et al. (2008) calculate a best fit of the surface to condensation temperature slope of 0.65 ^{∘}C ^{∘}C^{−1} in the ERA40 data, consistent with previous work that found a slope of 0.67 ^{∘}C ^{∘}C^{−1} (Connolley, 1996; Jouzel and Merlivat, 1984). The spread of condensation temperatures in the higherresolution MAR model suggests colder condensation temperatures than in the lowerresolution reanalysis (MassonDelmotte et al., 2008). The strength of the Antarctic inversion diminishes with increasing surface temperature (Connolley, 1996), and relatively warm Antarctic surface temperatures (e.g., $>\mathrm{20}$ ^{∘}C) are associated with condensation temperatures colder than the surface temperature (MassonDelmotte et al., 2008).
We analyze the relationship between surface temperature and condensation temperature in monthly MERRA2 reanalysis from 2008 through 2017 (Gelaro et al., 2017). We show the climatological zonalmean vertical profile of air temperature, the sum of the convective and largescale 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 takeaway is that the MERRA2 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 −10 ^{∘}C) 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.71–$\mathrm{0.75}\phantom{\rule{0.125em}{0ex}}\frac{{}^{\circ}\mathrm{C}}{{}^{\circ}\mathrm{C}}$ 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.
Using our nonlinear temperature reconstruction method, we model the condensation temperature for every pair of δ^{18}O and δD samples in the MD08 and GNIP datasets (MassonDelmotte et al., 2008; IAEA, 2001) 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 MERRA2 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 MERRA2 (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 bestfit 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 belowfreezing source evaporation is included (see below), in good agreement with previous Antarctic observational studies (Connolley, 1996; Jouzel and Merlivat, 1984). These slopes sit well within the range found in the MERRA2 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 bestfit 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 ${T}_{\mathrm{c}}=\mathrm{0.69}\frac{{}^{\circ}\mathrm{C}}{{}^{\circ}\mathrm{C}}\phantom{\rule{0.25em}{0ex}}{T}_{\mathrm{s}}\mathrm{8.2}\phantom{\rule{0.125em}{0ex}}{}^{\circ}\mathrm{C}$ as our base estimate to reconstruct Antarctic surface temperatures,; however, we consider an uncertainty of $\pm \mathrm{0.02}\frac{{}^{\circ}\mathrm{C}}{{}^{\circ}\mathrm{C}}$ in the slope. Our base estimate leads to good agreement with the observed relationship between global δ^{18}O and surface temperature (Fig. A10). We show our temperaturedependent isotope slopes in Fig. A11.
A3.3 Seasonality
Does seasonality in the hydrological cycle systematically bias climatological information in icecore waterisotope records? While the difference between precipitationweighted surface temperature and annual mean surface temperature is often discussed, this is not strictly the relevant comparison from the perspective of waterisotope ratios of precipitation. As discussed above, the critical comparison is between the annual mean surface temperature and the condensationweighted temperature, integrated over both time and altitude. Our analysis of the MERRA2 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 precipitationweighted 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 nonseasonal 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 ERAInterim reanalysis as well as global precipitation in the MERRA2 reanalysis. While the annual average Antarctic surface temperature and the precipitationweighted surface temperature are often different in either reanalysis product, we find little systematic bias. Across the Antarctic the monthtomonth and yeartoyear 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 precipitationweighted 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 precipitationweighted 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 condensationweighted 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 T_{0}? We examine the seasonality of Southern Hemisphere evaporation in the monthly MERRA2 reanalysis by comparing the annual average oversea surface temperatures to the evaporationweighted annual temperatures. Between 35 and 65^{∘} S, the bulk of Antarctic moisture sources and evaporationweighted 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.
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 Merlivat, 1984) as above, ${S}_{\mathrm{i}}=a+b\times 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 δ^{18}O in vapor and precipitation. We tune the model to yield the observed relationship between δD and δ^{18}O in global precipitation, rather than the relationship between δ values and environmental variables such as surface temperature.
Our target observational dataset includes waterisotope measurements of precipitation and surface snow from Antarctica and around the globe. The bulk of this compilation is that published by MassonDelmotte et al. (2008). We include additional published surface snow and precipitation measurements from the GNIP database (IAEA, 2001), surface traverses at Dome A (Xiao et al., 2013; Pang et al., 2015), Dahe et al. (1994), and a ^{17}O_{xs} 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 L2120i analyzer). Data are reported relative to the VSMOW (Vienna Standard Mean Ocean Water) standard and normalized to SLAP.
The global relationship between δD and δ^{18}O has an approximate slope of 8, as codified in the historical definition of the deuterium excess parameter (Dansgaard, 1964). However, the slope is not fundamental (Craig, 1961); 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 $\mathit{\delta}{{}^{\prime}}^{\mathrm{18}}\mathrm{O}$ in a global dataset for precipitation. They use a secondorder 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 δ^{18}O, given any amount of distillation, depends on the ratio of $\frac{\text{exp}{(}^{D}{\mathit{\alpha}}_{\mathrm{tot}})}{\text{exp}{(}^{\mathrm{18}}{\mathit{\alpha}}_{\mathrm{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 secondorder polynomial fit between $\mathit{\delta}{{}^{\prime}}^{\mathrm{18}}$O and δ^{′}D in our larger dataset compared to those found by Uemura et al. (2012): $A=\mathrm{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=\mathrm{28.5}$ and B=8.47) to define a logarithmic deuterium excess parameter, d_{ln}. We find no benefit or justification for using higherorder fits to this dataset.
We run SWIM to produce an ensemble of temperature trajectories representing a wide range of possible evaporation and condensation temperatures (T_{0} varies from 0 to 28 ^{∘}C; T_{c} from 27 to −60 ^{∘}C). We then compare the resulting cloud of modeled $\mathit{\delta}{{}^{\prime}}^{\mathrm{18}}\mathrm{O}$, δ^{′}D, and d_{ln}, finding a secondorder polynomial fit between the modeled δ^{′}D and $\mathit{\delta}{{}^{\prime}}^{\mathrm{18}}\mathrm{O}$ 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 $\mathit{\delta}{{}^{\prime}}^{\mathrm{18}}\mathrm{O}$ (Eq. 4). This is easily visualized in a plot of $\mathit{\delta}{{}^{\prime}}^{\mathrm{18}}\mathrm{O}$ and d_{ln} (e.g., Fig. A12a), as the average $\mathit{\delta}{{}^{\prime}}^{\mathrm{18}}\mathrm{O}$tod_{ln} relationship is flat in measured samples by definition.
Using the local closure assumption we find an optimal tuning of ${S}_{\mathrm{i}}=\mathrm{1}b\times T$ for b=0.00525 ^{∘}C^{−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 $\mathit{\delta}{{}^{\prime}}^{\mathrm{18}}\mathrm{O}$ 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 secondorder polynomial fit between modeled δ^{′}D and $\mathit{\delta}{{}^{\prime}}^{\mathrm{18}}\mathrm{O}$ to be identical to the observed fit: the observational target data represent a variety of timescales from subseasonal 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 higherelevation, colder Antarctic sites likely have more equatorward MSDs (Fig. 1), we should expect more depleted $\mathit{\delta}{{}^{\prime}}^{\mathrm{18}}\mathrm{O}$ in the target data to be associated with slightly warmer T_{0} (that is, modeled results from a single value of T_{0} should not be strictly flat in the $\mathit{\delta}{{}^{\prime}}^{\mathrm{18}}\mathrm{O}$–d_{ln} 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 T_{0} between $\mathit{\delta}{{}^{\prime}}^{\mathrm{18}}\mathrm{O}$ values of −40 ‰ and −55 ‰ in the observational dataset. This is a fairly small shift in mean T_{0} compared to the range of T_{0} that may contribute to a site but should give some downward curvature to model results of equal T_{0} at the coldest T_{c} values.
The appropriate weighting of model realizations with different T_{0} 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 T_{0}, 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., $\mathrm{0.005}\le b\le \mathrm{0.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.
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, temperaturedependent 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, ${S}_{\mathrm{i}}=ab\times Tc\times {T}^{\mathrm{2}}$. 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 S_{i} with decreasing temperature could be plausible. Only very small values of the secondorder term, c, those of order $\mathrm{5}\times {\mathrm{10}}^{\mathrm{6}}$ ^{∘}C^{−2}, are reconcilable with the observed $\mathit{\delta}{{}^{\prime}}^{\mathrm{18}}\mathrm{O}$tod_{ln} 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 moistureweighted average of a set of independent pseudoadiabatic 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 moistureweighted average of precipitation from unmixed pathways. There are three processes associated with mixing to consider.

Nonuniqueness. Parcels that experienced different evaporation conditions can arrive at a condensation site of a given temperature with different isotopic values.

Mixinginduced 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 moistureweighted average of unmixed pathways.

Nonlinear mixing. For isobaric mixing of equalmassed air parcels, the final mixed temperature reflects their massweighted 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 moistureweighted average of the unmixed parcels.
Moistureweighted 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 (massweighted) 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 ERAInterim reanalysis (Dee et al., 2011) over the Southern Hemisphere oceans. In summer, the daytoday 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 synopticscale variability such as lagged 2 or 5 d temperature differences or gridpointto gridpoint 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 longterm 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 (T_{0}=10 ^{∘}C), 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 −30 ^{∘}C. 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 highend estimate based on the above analysis of daily temperature differences). This process is repeated 10 000 times to create a distribution of final δ^{18}O (Fig. A14a), δD (Fig. A14b), and d_{ln} (Fig. A14c) values of precipitation at −30 ^{∘}C, 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 moistureweighted 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 T_{0} and T_{c}, though the differences in moistureweighted 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 moistureweighted mean isotopic value. In Fig. A14e–h we show the distribution of the isotopic value of precipitation at −30 ^{∘}C 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 moistureweighted mean of the unmixed parcels. Differences between the mean isotopic values of mixed and unmixed parcels are less than 0.2 ‰ for δ^{18}O, 1.5 ‰ for δD, and 0.01 ‰ for d_{ln}. 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 moistureweighted 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 δ^{18}O and condensation temperature and the relationship between d_{ln} and evaporation temperature are similar whether or not mixing is present. Below we investigate the influence of mixing on our temperature reconstruction technique.
A6 Optimal coordinates for reconstruction technique
Consider a water sample with mean values of δ^{18}O 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 d_{xs} and d_{ln} 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 T_{0} using coordinates of δ^{18}O and δD (Fig. A15a), δ^{18}O and d_{xs} (Fig. A15b), and δ^{18}O and d_{ln} (Fig. A15c). The uncertainties in the position of the measurement along both the x axis (δ^{18}O) and y axis (either δD, d_{xs}, or d_{ln}) combine to give the total uncertainty in the position on the T_{0} surface, shown as the targets in Fig. A15a–c. The total combined uncertainty in the estimation of T_{0} is shown as probability density functions (PDFs) for each method in Fig. A15d. All estimates yield the same mean value of reconstructed T_{0}; 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 d_{xs} and d_{ln} reconstructions are similar, the d_{ln} reconstruction has a narrower PDF and thus more confident reconstruction. This is because the T_{0} isotherms are most separated along and most perpendicular to the d_{ln} axis. The advantage of the separation of isotherms along the d_{ln} 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 d_{xs} at more depleted δ^{18}O values, that uncertainty will be spread across a wider range of isotherms. The axis of the influence of T_{0} on d_{xs} (and δD) is rotated with respect to its axis of variability in d_{xs}.
This result is of course ultimately tied to the same reasons that d_{ln} provides a better qualitative proxy for source region changes than d_{xs}. The initial imprint of the source conditions is better preserved in d_{ln} than d_{xs}. 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 δ^{18}O and d_{ln} measurements (corrected for changes in seawater δ^{18}O) from the eight Antarctic deep icecore records examined in this study.
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 200year 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 sampletosample variability is larger than the relative uncertainty, making the shading difficult to see at the scale plotted here.
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.
We compare the nonlinear temperature reconstructions of all eight icecore 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 T_{0} and T_{s} using the linear and nonlinear technique (Fig. A21). Reconstructed site surface temperatures, T_{s}, 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, T_{0}, there is dramatic improvement in coherence amongst the records when using the nonlinear technique. The increase in shared variance (R^{2}) 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 waterisotope 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 d_{ln} is a more faithful qualitative proxy for source region conditions than the linearly define d_{xs}. Compare the correlation matrices of the excess parameters in Fig. A23. The nonlinearly reconstructed T_{0} and d_{ln} parameters share the same correlation pattern amongst the ice cores and show substantially more coherence than either linearly reconstructed T_{0} or d_{xs}. The correlation pattern of d_{ln} and d_{xs} between all core sites (Fig. A24) reveals how nonlinear effects alter the traditionally defined d_{xs} at the coldest Antarctic temperatures. The broad change from positive to negative correlation of d_{ln} to d_{xs} across sites is a reflection of the change in sign of $\frac{\partial {d}_{\mathrm{xs}}}{\partial {T}_{\mathrm{c}}}$ as a function of T_{c}.
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 T_{c} and T_{0} from a set of δ^{18}O and d_{ln} 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 T_{c} and T_{0} 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 δ^{18}O and d_{ln}. Thus, uncertainty arising from a given parameter may vary between icecore sites. In general, varying a parameter in the model results in patterns of the partial slopes in δ^{18}O and d_{ln} with T_{c} and T_{0} that are similar, but shifted in the T_{c} and T_{0} space. A consequence of this is that the uncertainty in absolute values of reconstructed T_{c} and T_{0} is generally larger than uncertainty in their relative variability.
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.0054 ^{∘}C^{−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 ${\mathit{\alpha}}_{\mathrm{diff}}^{\mathrm{18}}$ 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 endmember 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 T_{0} up to 1.5 ^{∘}C for the WDC record, while differences in absolute T_{c} are smaller (≤1 ^{∘}C). Relative variability in T_{0} and T_{c} is similar when using either closure assumption and ≤0.3 ^{∘}C.
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 ERAInterim reanalysis. We can also ask how different our reconstructed T_{c} and T_{0} in WDC would be if we used fixed mean values of initial relative humidity rather than values that depend on T_{0}. 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 T_{0} and T_{c} as functions of δ^{18}O and d_{ln} 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 −50 ^{∘}C. We generate random pairs of pseudoadiabatic cooling pathways ending at every value of T_{c} and random values of T_{0} 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 T_{c}. To mix, parcel temperatures must be above an absolute threshold temperature (−15 ^{∘}C) 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 T_{c} between −20 and −50 ^{∘}C in increments of 0.1 ^{∘}C (a total of 1.5×10^{4} unmixed and 7.5×10^{5} mixed cooling paths). We then interpolate the results of both the mixed pathways and the moistureweighted averages of the unmixed pathways to create maps of T_{0} and T_{c} as functions of δ^{18}O and d_{ln} (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 T_{c} and T_{0} maps resulting from the mixed (denoted with the subscript M) and unmixed pathways (Fig. A26). In Fig. A26 we show the histograms of $\mathrm{\Delta}{T}_{\mathrm{c}}={T}_{\mathrm{c},\mathrm{M}}{T}_{\mathrm{c}}$ and $\mathrm{\Delta}{T}_{\mathrm{0}}={T}_{\mathrm{0},\mathrm{M}}{T}_{\mathrm{0}}$. We test a range of mixing and threshold values in multiple Monte Carlo simulations. In all cases the mean values of ${T}_{\mathrm{c},\mathrm{M}}{T}_{\mathrm{c}}$ and ${T}_{\mathrm{0},\mathrm{M}}{T}_{\mathrm{0}}$ are very near zero ($<\pm \mathrm{0.06}$ ^{∘}C for T_{c}, and $<\pm \mathrm{0.02}$ ^{∘}C for T_{0}). 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.2 ^{∘}C for T_{c} and ±0.35 ^{∘}C for T_{0} in all tests. We find similar results when using isobaric rather than pseudoadiabatic cooling pathways. Including moisture sources with ${T}_{\mathrm{0}}<\mathrm{0}{}^{\circ}$C in our mixing analysis has no significant influence on the mean difference between the weighted mean maps of either T_{c} or T_{0}, though expanding the range of moisture sources to include T_{0}<0 ^{∘}C does increase the range of ${T}_{\mathrm{0},\mathrm{M}}{T}_{\mathrm{0}}$ by over a degree, in agreement with the analysis of unmixed pathways.
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 T_{c} and T_{0} reconstruction interpolated at each pair of δ^{18}O and d_{ln} 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 ${}^{\mathrm{18}}{a}_{\mathrm{diff}}=\mathrm{1.009}\pm \mathrm{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 T_{0} and RH_{0} as the mean absolute difference in reconstructions using climatological fits from the NCEP/NCAR and ERAInterim 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 T_{c} and T_{0} differences are nonuniform and unevenly populated, making interpolation in the δ^{18}O–d_{ln} 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 δ^{18}O and d_{ln} symmetrically around the base scenario. Finally, we include the additional uncertainty in our estimates of relative T_{s} variability arising from the T_{c}toT_{s} relationship outlined in Sect. A3.2. We use the mean absolute difference of reconstructions using ${T}_{\mathrm{c}}\propto \mathrm{0.69}\pm \mathrm{0.02}\frac{{}^{\circ}\mathrm{C}}{{}^{\circ}\mathrm{C}}\phantom{\rule{0.25em}{0ex}}{T}_{\mathrm{s}}$ to estimate this uncertainty, which is added in quadrature to the above uncertainties in T_{s}. An example of this spread of uncertainty in the WDC T_{s} reconstruction is shown in Fig. A27. Reconstructions of T_{s}, T_{c}, and T_{0} for several major icecore records along with the combined relative uncertainty are those reconstructions is shown in Fig. 9.
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 watertagged 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 δ^{18}O for every value of condensation temperature owing to the influence of evaporation temperature, there are not necessarily unique pairs of δ^{18}O and d_{ln} for every pair of T_{c} and T_{0} if T_{0} can be both above and below freezing. The T_{c} and T_{0} surfaces fold over on themselves in the δ^{18}O and d_{ln} space for values of T_{0} 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 nonuniqueness in the waterisotope state spaces. We run SWIM through a large field of T_{0} values from −28 to 28 ^{∘}C. We can in principle resolve the nonuniqueness problem by combining reconstructions from either side of the folded surface based on the contribution of total moisture represented by each pair of nonunique paths. Knowing that belowzero moisture sources contribute far less to the total moisture reaching Antarctic sites (Fig. A28), we examine two reasonable methods of moistureweighting the reconstructions. In the first we simply weight each pair of reconstructed temperatures by the final mixing ratio (r_{s(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 higherelevation, colder sites compared to our GCMbased MSD estimates. In this r_{s(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 % T_{0}>0 ^{∘}C, 10 % T_{0}<0 ^{∘}C). 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 nonuniqueness leads to very small differences in reconstructed T_{c} and T_{0} variability: the standard deviation of differences in reconstructed T_{c} is less than 0.07 ^{∘}C using either method and less than 0.19 ^{∘}C for T_{0}. Attempting to account for this nonuniqueness does, however, lead to persistent mean offsets in absolute temperature; in particular, we find colder absolute values of reconstructed T_{0} for all icecore sites.
While these results are interesting, this attempt to account for nonuniqueness likely does not actually improve the absolute temperature reconstructions. Given the shape of the folded temperature surfaces in the modeled δ^{18}O and d_{ln} space, as well as the actual values of δ^{18}O and d_{ln} in icecore measurements, the model must extrapolate to extremely cold T_{0} values for the below 0 ^{∘}C side of the folded surface. These values of T_{0} are far colder (>10 ^{∘}C 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. Nearsurface 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 T_{0} when a much warmer one (and perhaps a reduced RH_{0}) is actually correct. The net result of these considerations is that the analysis above should represent a quite conservative estimate of the influence of nonuniqueness 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 icecore records and compare to previously published linear reconstructions. We use records of δ^{18}O and δD measurements (and calculate d_{xs} and d_{ln}) 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 icecore δD and d_{xs} for the linear reconstruction and δ^{18}O and d_{ln} 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, T_{source}, equivalent to our evaporation temperature, T_{0}, and for the site surface temperature, T_{site}. We convert our reconstructed condensation temperatures, T_{c}, 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 waterisotope 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 waterisotope to temperature relationships is the deuterium excess parameters, $\frac{\partial {d}_{\mathrm{xs}}}{\partial {T}_{\mathrm{c}}}$ and $\frac{\partial {d}_{\mathrm{xs}}}{\partial {T}_{\mathrm{0}}}$. If one assumes these slopes are linear over a given range in T_{0} and T_{c} when in reality they are nonlinear, one will attribute a given change in Δd_{xs} 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 δ^{18}O 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 ±2 ^{∘}C (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 −5 ^{∘}C. 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 −5 ^{∘}C. 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 d_{xs} versus d_{ln} (Markle et al., 2017; Uemura et al., 2012) (see Sect. 1.2) and ultimately a consequence of the same distillation effects.
A11 Threeparameter reconstructions
In the approach outlined above, we have considered the boundary conditions T_{c} and T_{0} to be the only independent input variables. In particular, we have assumed that the source region relative humidity, RH_{0}, is a dependent variable whose value is not fixed but determined by climatological correlation with T_{0}. 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 RH_{0} is dependent on T_{0} and reconstruct three independent climate variables (T_{c}, T_{0}, and RH_{0}) if we have three independent constraints. While δ^{18}O and d_{ln} alone are not sufficient, ${}^{\mathrm{17}}{\mathrm{O}}_{\mathrm{xs}}=\mathit{\delta}{{}^{\prime}}^{\mathrm{17}}\mathrm{O}\mathrm{0.528}\mathit{\delta}{{}^{\prime}}^{\mathrm{18}}$O (Landais et al., 2008) can in principle provide the necessary additional information. We can allow T_{c}, T_{0}, and RH_{0}, to all vary as independent variables, defining a threedimensional parameter space through which SWIM is run to produce threedimensional isotope state spaces.
While promising, this method currently has practical limitations. Our model does not reproduce the observed ^{17}O_{xs}toδ^{18}O 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 ^{17}O_{xs} in Antarctic precipitation may be too large to offer useful discrimination amongst variations in T_{c}, T_{0}, and RH_{0} (Schoenemann et al., 2014).
The ^{17}O_{xs} 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 ^{17}O_{xs} lead to large changes in the absolute value of reconstructed source region conditions (T_{0} and RH_{0}). It is worth noting that absolute values of ^{17}O_{xs} are 3 orders of magnitude smaller than the deuterium excess. Further, preliminary testing suggests that there may be significant nonuniqueness to address; that is, a position in the threeparameter space, defined by δ^{18}O, d_{ln}, and ^{17}O_{xs}, 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), waterisotope 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.
The temperature reconstructions, underlying data, and Simple Water Isotope Model code are available through a public GitHub repository (doi: https://doi.org/10.5281/zenodo.6510097, Markle and Steig, 2022; https://github.com/bradleymarkle/simple_water_isotope_model, Markle, 2022) and the United States Antarctic Research Program (USAP) Data Center. The icecore data are all already available through links provided in the primary papers cited for each dataset.
BRM conceived of the study and, with EJS, designed the model and analyses. BRM and EJS wrote the paper.
The contact author has declared that neither they nor their coauthor 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 SteenLarsen, 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.
This research has been supported by the National Science Foundation (grant nos. 1043092, 1141839, and 1443105).
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, https://doi.org/10.1029/2019GL082965, 2019. a, b, c, d
Bakhshaii, A. and Stull, R.: Saturated PseudoadiabatsA Noniterative Approximation, J. Appl. Meteorol. Climatol., 52, 5–15, 2013. a
Barkan, E. and Luz, B.: High precision measurements of ${}^{\mathrm{17}}\mathrm{O}{/}^{\mathrm{16}}\mathrm{O}$ and ${}^{\mathrm{18}}\mathrm{O}{/}^{\mathrm{16}}\mathrm{O}$ ratios in H_{2}O, Rapid Commun. Mass Spect., 19, 3737–3742, 2005. a
Barkan, E. and Luz, B.: Diffusivity fractionations of H${}_{\mathrm{2}}^{\mathrm{16}}$O/H ${}_{\mathrm{2}}^{\mathrm{17}}$O and H${}_{\mathrm{2}}^{\mathrm{16}}$O/H${}_{\mathrm{2}}^{\mathrm{18}}$O in air and their implications for isotope hydrology, Rapid Commun. Mass Spect., 21, 2999–3005, 2007. a, b
Benetti, M., SteenLarsen, 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 icesheet dynamics and the onset of 100,000year 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 MERRA2 reanalysis, J. Climate, 30, 1177–1196, 2017. a
Braconnot, P., Harrison, S. P., Kageyama, M., Bartlein, P. J., MassonDelmotte, V., AbeOuchi, A., OttoBliesner, 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 millennialscale 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., GotoAzuma, K., Kawamura, K., Fujita, S., Motoyama, H., Hirabayashi, M., Uemura, R., Stenni, B., Parrenin, F., He, F., Fudge, T. J., and Steig, E.: Abrupt iceage 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., SherriffTadano, S., Ritz, C., Lefebvre, E., Edwards, J., Kawamura, K., Oyabu, I., Motoyama, H., Kahle, E. C., Jones, T. R., AbeOuchi, 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., OttoBliesner, 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 TransAntarctica 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, https://doi.org/10.3402/tellusa.v16i4.8993, 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., MongeSanz, B. M., Morcrette, J.J., Park, B.K., Peubey, C., de Rosnay, P., Tavolato, C., Thépaut, J.N., and Vitart, F.: The ERAInterim reanalysis: Configuration and performance of the data assimilation system, Q. J. Roy. Meteorol. Soc., 137, 553–597, 2011. a, b
EPICA Community Members: Onetoone coupling of glacial climate variability in Greenland and Antarctica, Nature, 444, 195–198, https://doi.org/10.1038/nature05301, 2006. a
Epstein, S., Sharp, R. P., and Goddard, I.: Oxygenisotope 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 modernera retrospective analysis for research and applications, version 2 (MERRA2), 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.: Firstprinciples diffusivity ratios for kinetic isotope fractionation of water in air, Geophys. Res. Lett., 47, e2020GL089999, https://doi.org/10.1029/2020GL089999, 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, https://doi.org/10.1029/2009JD012384, 2010. a, b, c, d, e
IAEA: GNIP Maps and Animations, International Atomic Energy Agency, Vienna, http://isohis.iaea.org (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, https://doi.org/10.1029/2020JD033300, 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 40year 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 δ^{18}O and δD in Antarctic precipitation revisited, Global Biogeochem. Cy., 17, 1017, https://doi.org/10.1029/2002GB001910, 2003. a, b, c, d, e, f, g, h
Kurita, N., Hirasawa, N., Koga, S., Matsushita, J., SteenLarsen, H. C., MassonDelmotte, V., and Fujiyoshi, Y.: Influence of largescale 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, https://doi.org/10.1073/pnas.1618374114, 2017. a
Landais, A., Barkan, E., and Luz, B.: Record of δ^{18}O and ^{17}Oexcess in ice from Vostok Antarctica during the last 150,000 years, Geophys. Res. Lett., 35, L02709, https://doi.org/10.1029/2007GL032096 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, https://doi.org/10.1038/NGEO411, 2009. a
Markle, B.: Simple Water Isotope Model, GitHub [code], https://github.com/bradleymarkle/simple_water_isotope_model, last access: 3 May 2022. a
Markle, B. R. and Steig, E. J.: bradleymarkle/simple_water_isotope_model: Simple Water Isotope Model (Version v0), Zenodo [data set], https://doi.org/10.5281/zenodo.6510097, 2022. a
Markle, B., Bertler, N., Sinclair, K., and Sneed, S.: Synoptic variability in the Ross Sea region, Antarctica, as seen from backtrajectory modeling and ice core analysis, J. Geophys. Res.Atmos., 117, D02113, https://doi.org/10.1029/2011JD016437, 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 DansgaardOeschger 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 highlatitude aerosols, water isotopes and the hydrologic cycle, Nat. Geosci., 11, 853–859, https://doi.org/10.1038/s4156101802109, 2018. a, b
MassonDelmotte, V., Kageyama, M., Braconnot, P., Charbit, S., Krinner, G., Ritz, C., Guilyardi, E., Jouzel, J., AbeOuchi, A., Crucifix, M., Gladstone, R. M., Hewitt, C. D., Kitoh, A., LeGrande, A. N., Marti, O., Merkel, U., Motoi, T., Ohgaito, R., OttoBliesner, 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 icecore constraints, Clim. Dynam., 26, 513–529, 2006. a, b
MassonDelmotte, 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 H${}_{\mathrm{2}}^{\mathrm{16}}$O, HD^{16}O, and H${}_{\mathrm{2}}^{\mathrm{18}}$O in gases, The J. Chem. Phys., 69, 2864–2871, 1978. a, b, c
Merlivat, L. and Jouzel, J.: Global climatic interpretation of the deuteriumoxygen 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 solidevapeur et liquidevapeur 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, https://doi.org/10.1175/15200450(1967)006<0203:OTCOSV>2.0.CO;2, 1966. a
OttoBliesner, 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., MassonDelmotte, V., Prie, F., SteenLarsen, H. C., Risi, C., Li, Y., Jouzel, J., Wang, Y., He, J., Minster, B., and Falourd, S.: Spatial distribution of ^{17}Oexcess 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 CraigGordon model, J. Geophys. Res.Atmos., 114, D20108, https://doi.org/10.1029/2009JD012054, 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., MassonDelmotte, V., and Vimeux, F.: Understanding the ^{17}O excess glacialinterglacial variations in Vostok precipitation, J. Geophys. Res.Atmos., 115, D10112, https://doi.org/10.1029/2008JD011535, 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 9780549142577, 2007. a
Schmidt, G., Bigg, G., and Rohling, E.: Global seawater oxygen18 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 waterisotopologue record from WAIS Divide, Antarctica: Controls on glacialinterglacial changes in ^{17}Oexcess 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, https://doi.org/10.1029/2009GL040242, 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 icesheet 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.: Continuousflow analysis of δ^{17}O, δ^{18}O, and δD of H_{2}O on an ice core from the South Pole, Front. Earth Sci., 9, 640292, https://doi.org/10.3389/feart.2021.640292, 2021. a
Stenni, B., Jouzel, J., MassonDelmotte, 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 lateglacial highresolution 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., MassonDelmotte, 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., LemieuxDudon, B., Maggi, V., MassonDelmotte, 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 seesaw 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, https://doi.org/10.1029/2008JD010209, 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, https://doi.org/10.1029/2009GL041960, 2010. a, b, c, d, e, f, g
Uemura, R., MassonDelmotte, V., Jouzel, J., Landais, A., Motoyama, H., and Stenni, B.: Ranges of moisturesource temperature estimated from Antarctic ice cores stable isotope records over glacial–interglacial cycles, Clim. Past, 8, 1109–1125, https://doi.org/10.5194/cp811092012, 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, https://doi.org/10.1016/S0012821X(02)009500, 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, https://doi.org/10.1038/nature12376, 2013. a, b, c
WAIS Divide Project Members: Precise interpolar phasing of abrupt climate change during the last ice age, Nature, 520, 661–665, https://doi.org/10.1038/nature14401, 2015. a
Werner, M., Jouzel, J., MassonDelmotte, 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., MassonDelmotte, V., Stenni, B., and Jouzel, J.: Deglaciation records of ^{17}Oexcess in East Antarctica: reliable reconstruction of oceanic normalized relative humidity from coastal sites, Clim. Past, 8, 1–16, https://doi.org/10.5194/cp812012, 2012. a
Xiao, C., Ding, M., MassonDelmotte, 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 SST_{0} and RH_{0} for a given T_{0} determined by a spline fit to NCEP/NCAR climatology, and a tuned supersaturation such that $S=\mathrm{1}\mathrm{0.00525}T$, where T is in degrees Celsius (^{∘}C). See the Appendix for details.
 Abstract
 Introduction
 A (relatively) simple waterisotope model
 Temperaturedependent slopes
 Nonlinear temperature reconstructions
 Discussion
 Conclusions
 Appendix A: Simple Water Isotope Model
 Code and data availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References
 Abstract
 Introduction
 A (relatively) simple waterisotope model
 Temperaturedependent slopes
 Nonlinear temperature reconstructions
 Discussion
 Conclusions
 Appendix A: Simple Water Isotope Model
 Code and data availability
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References