Improving temperature reconstructions from ice-core water-isotope records

. 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 origi-nates. Temperature reconstructions obtained from ice-core records generally rely on linear approximations of the relationships among local temperature, source temperature, and water-isotope values. However, there are important nonlinearities that signiﬁcantly affect such reconstructions, particularly for source region temperatures. Here, we describe a relatively simple water-isotope distillation model and a novel temperature reconstruction method that accounts for these nonlinearities. Further, we examine in detail many of the parameters, assumptions, and uncertainties that underlie water-isotope distillation models and their inﬂuence on these temperature reconstructions. We provide new reconstructions of absolute surface temperature, condensation temperature, and source region evaporation temperature for all long Antarctic ice-core records for which the necessary data are available. These reconstructions differ from previous estimates due to both our new model and reconstruction technique, the inﬂuence of which is investigated directly. We also provide thor-ough uncertainty estimates for all temperature histories. Our reconstructions constrain the pattern and magnitude of polar ampliﬁcation in the past and reveal asymmetries in the temperature histories of East and West Antarctica. . 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.


Introduction
Stable-isotope ratios of water have been the foundational proxy for polar paleoclimate research for more than a half-century (Langway, 1958;Gonfiantini, 1959;Dansgaard, 1964). Primarily used as a temperature proxy, stratigraphic records of water-isotope ratios in ice sheets provide detailed histories of Earth's climate over hundreds of thousands of years (Dansgaard et al., 1969;Petit et al., 1999), providing insight into the past magnitudes, spatial patterns, and phasing of climate change across the globe EPICA Community Members, 2006;WAIS Divide Project Members, 2013, 2015. Both oxygen and hydrogen have stable isotopes whose ratios ( 18 O/ 16 O and 2 H/ 1 H) are commonly expressed as deviations, δ 18 O and δD, from Vienna Standard Mean Ocean Water (VSMOW) in per mill (‰): 1322 B. R. Markle and E. J. Steig: Improving temperature reconstructions from ice cores 1999; Jouzel et al., 1997;Masson-Delmotte et al., 2008). Air temperatures during condensation (Petit et al., 1999;Jouzel et al., 1997) and during initial moisture evaporation (Vimeux et al., 2002) can be reconstructed from ice-core water-isotope records if the relevant scaling relationships can be determined from theory, models, or observations (Vimeux et al., 2002;Kavanaugh and Cuffey, 2002;Stenni et al., 2010).
Here, we examine the widely used assumption of linearity in the scaling relationships between water-isotope ratios and temperature.

Temperature reconstructions
Any interpretation of water-isotope ratios as a proxy for temperature requires a model, whether conceptual, statistical, or numerical. A conceptual model of progressive distillation and integrated fractionation (e.g., Dansgaard, 1964) is sufficient to qualitatively interpret variations in water-isotope ratios as variations in temperature in the high latitudes. The simplest quantitative interpretation of ice-core water-isotope records relies on the empirical correlation between observed water-isotope ratios of precipitation and surface temperature at the precipitation site (Petit et al., 1999;Jouzel et al., 1997). A limitation of this approach is the possibility to conflate the "spatial slope" between water isotopes and temperature, which is the relationship observed across a range of modern sites, and the "temporal slope", which is the relationship at a single point through time (Jouzel et al., 1997). Relatedly, this approach also does not account for simultaneous and independent changes in evaporation conditions, which can impact high-latitude water-isotope ratios in several ways. Initial evaporation temperature, together with the condensation temperature, determines the total temperature gradient through which moisture must be distilled to reach a given site. Further, evaporative conditions set the initial isotopic values of the vapor before distillation. The isotope ratios of vapor above the ocean depend on the temperature during evaporation, the isotopic values of the seawater, and the occurrence of kinetic fractionation during evaporation, which is driven by sub-equilibrium relative humidity and influenced by sea surface temperature and wind speed (Merlivat and Jouzel, 1979;Jouzel et al., 1982). A more complete approach to reconstructing temperature from water-isotope records is to employ numerical models that account for the combined influence of variability in both evaporation and condensation temperatures, as well as other factors. Reconstructing two unknowns (i.e., both evaporation source and condensation site temperatures) requires two constraints, which are provided by the oxygen and hydrogen isotope ratios and the relationship between them. While the oxygen and hydrogen isotope systems have similar behavior in the atmosphere, there are differences in their response to the same environmental conditions and to processes such as kinetic fractionation. The deuterium excess is the weighted difference between δ 18 O and δD, d xs = δD−8 × δ 18 O, and is commonly used to quantify these differences (Dansgaard, 1964;Merlivat and Jouzel, 1979).
Changes in water-isotope parameters measured in precipitation at an ice-core 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 water-isotope distillation models for the ice-core 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).

Nonlinearities in isotope fractionation and the deuterium excess definition
The temperature reconstruction approach described above depends on the assumption that the parameters, β and γ , are fixed in time and independent of temperature. However, the β and γ parameters, as diagnosed from model simulations, are found to be different for different ice-core sites, whose only distinguishing characteristics are differing modern surface conditions (e.g., Stenni et al., 2010;Uemura et al., 2012). Thus, β and γ depend on the site conditions, which obviously change over time.
Another issue with the linear reconstruction approach is the definition of the deuterium excess parameter (Uemura et al., 2012;Markle et al., 2017). The origin of the slope in the definition of deuterium excess is an empirical fit to global precipitation measurements (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 second-order polynomial to a compilation of δ 18 O and δ D data, where δ x = ln (1 + δ x ), and defined a phenomenological, nonlinear deuterium excess parameter: with coefficients A = −28.5 and B = 8.47 (note that the coefficients and δ values are unitless; for example, δ 18 O = 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 water-isotope-based temperature reconstructions and suggest an improved technique.

A (relatively) simple water-isotope model
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 high-latitude water-isotope ratios are a useful temperature proxy. They are driven by two basic nonlinear processes, the Clausius-Clapeyron relationship and Rayleigh distillation (see Sect. A2.2).
Other processes can be important as well. The temperature dependence of fractionation factors, for example, generally amplifies the temperature relationship. While any single precipitation event at a site may be subject to a variety of additional factors and processes, the long-term mean is strongly influenced by climatological moisture distillation. 1 Our model distills moisture down thermodynamic pathways defined by temperature and pressure. Changes in 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 pseudo-adiabatic, consistent with isentropic moisture transport to the Antarctic (Bailey et al., 2019) and the basic assumption of Raleigh distillation that moisture is removed after precipitation. A superposition of many thermodynamic pathways is required to represent a single Antarctic precipitation site, reflecting both the range of precipitation conditions experienced at a site and moisture transport from sources with a distribution of evaporative conditions (Markle et al., 2017, Fig. 1). An example of a set of these pathways is shown in Fig. 2. We use climatological correlations to relate initial evaporation air temperature, 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 non-uniqueness. We also investigate the influence of mixing during both evaporation and transport, the potential influence of seasonality (and intermittency) of precipitation and evaporation, and the relationship between surface temperature at a site and the vertically integrated, precipitation-weighted condensation temperature above that 1 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.  (Markle et al., 2017). (a) Difference, in degrees latitude, between the latitude of precipitation and mean latitude of evaporation (effectively mean transport in degrees of latitude to any site). (b) Mean latitude of evaporation vs. latitude of precipitation. All longitudes are shown by the black line; longitudes encompassing West Antarctica are shown by the blue line, while longitudes encompassing East Antarctica are shown in red. (c) The mean evaporation latitude (in • S) of precipitation falling at all Antarctic grid points. Select ice-core sites shown in white. (d) The relationship between mean evaporation latitude and site elevation across Antarctica. Select ice-core locations shown in color; see text for site details. (e) Latitudinal moisture source distributions for select Antarctic ice-core sites, colored by site elevation. site. We make no attempt to model post-depositional processes.

Temperature-dependent slopes
To investigate the sensitivity of water-isotope ratios of Antarctic precipitation to site and source conditions, we use SWIM to model isotopic state spaces. We run SWIM through a large ensemble of temperature pathways defined by T = T 0 , T 0 − dT , . . ., T c , with dT = 0.1 • C. We run nearly 24 000 trajectories to fill the plausible parameter space of 0 • C ≤ T 0 ≤ 28 • C and −70 • C ≤ T c ≤ 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 2 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 = 1 − 0.00525T , where T is in degrees Celsius ( • C). See the Appendix for details. 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)  and ∂d xs ∂T 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 ∂d xs 1326 B. R. Markle and E. J. Steig: Improving temperature reconstructions from ice cores 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 ∂d xs ∂T c (Fig. 4c) can actually be seen directly in Antarctic ice-core 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 ∂d xs ∂T 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 ice-core sites. Using the same isotope model, they calculated different scalings for each site. However, by assuming that these slopes are constant for each site, they do not consider the possibility that one site's conditions may have been more like another's in the past. Recognizing this as well as the inability of their model to simultaneously match observed site temperature and δ values, Uemura et al. (2012) create several reconstructions for the Dome Fuji site utilizing different linearizations of the model. They do not, however, attempt a reconstruction that accounts for the nonlinearities in the water-isotope-temperature relationships.
The solution to this issue of slope nonlinearity, within the linear isotope temperature reconstruction framework (Eqs. 2 and 3), is not obvious. The nonlinearities in the slopes of the isotope surfaces depend on the absolute condensation and evaporation temperatures that represent the target of the reconstruction, which are of course not known a priori. We next present a novel temperature reconstruction framework, which takes into account the inherent nonlinearities in the water-isotope fractionation process.

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 water-isotope model around the modern climate state, this method accounts for the changes in the slopes that depend on the mean state. Further, there is no need to find analytical solutions to the model or fit families of high-order polynomials to the results.
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 second-order 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.

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 tem-perature. 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 MERRA-2 reanalysis data (Gelaro et al., 2017). We also examine the seasonality and vertical distribution of Antarctic precipitation and reevaporation. Based on these analyses we use a simple, linear temperature-dependent relationship to estimate weighted, annual mean surface air temperature from our condensation temperature reconstructions and account for the uncertainty in this relationship.
We examine in detail how assumptions and modeling choices concerning initial evaporation (Appendix A2.1) impact our results, as well as the potential for bias arising from seasonality in evaporation from the ocean. We also examine the influence of non-uniqueness on our temperature reconstruction technique arising from below-freezing evaporation conditions (Appendix A9.4).
Finally, we conduct an extensive uncertainty analysis of our temperature reconstructions (Appendix A9). We calculate numerous isotope state spaces from the same temperature parameter space using multiple iterations of the model in which model assumptions are altered and parameters are varied over plausible ranges. We calculate the uncertainty in In the time series plots, each record is colored by its Holocene average δ 18 O value. See Sect. 4.2 in the text for details and references for the ice-core records. All original records are interpolated to even 50-year time spacing on the Buizert et al. (2018) synchronized timescale where possible, or they are plotted on original published timescales. All records are ordered by their approximate modern surface temperature in the cross-correlation matrices. 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 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 water-isotope ratios recorded in the ice cores, leaving the identification of the causes of those changes in temperature for future work. This means that neither the sites of precipitation nor evaporation are strictly fixed in space. Indeed, variability in moisture source locations may cause a significant amount of variation in reconstructed evaporation temperature (Markle et al., 2017). We aim to maintain the broad utility of our reconstructions by building in as few assumptions about the sources of the temperature variability as possible.

Comparison of linear and nonlinear reconstruction techniques
Linear temperature reconstruction using SWIM We evaluate the significance of our approach by comparing our nonlinear reconstructions of condensation and evaporation temperature with reconstructions following the traditional linear approach. We calculate the equivalent linear β and γ coefficients for each of eight ice-core sites by regressing the 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 ice-core location using the two sets of linear β and γ coefficients found above and the traditional linear method. We then compare these linear reconstructions to our full nonlinear reconstruction and show the difference in reconstructed surface temperature and evapora- tion temperature in Fig. 10. We also show the residuals between the two linear reconstructions. In general, the residuals are not constant offsets but vary as a function of reconstructed temperature, demonstrating the temperature dependence of the slopes. They also differ, sometimes in sign, between the ice-core sites (Fig. A20). Linearization can obscure true variability or introduce spurious variability into the reconstructions, depending on the actual conditions of the site over time.
In the case of the surface temperature reconstruction, the errors introduced by linearization can be up to ±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 ice-core site surface temperature. The errors introduced into the reconstructed evaporation temperatures by ignoring the nonlinearities can be nearly as large as the total reconstructed variability. It is thus problematic to attempt reconstructing evaporation temperatures without accounting for nonlinearities.
In Appendix A8 we show that using the nonlinear reconstruction technique yields greater correlation amongst all records of 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 ice-core site temperature reconstructions from the linear and nonlinear reconstruction techniques have relatively small differences. However, as seen in Fig. 10, there are small artifacts arising from slight nonlinearity in the δ 18 O-to-T 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 liquid-only and ice-only 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 O-to-T 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 water-isotope models (Stenni et al., 2004Uemura 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. Reconstructions of (b) moisture source evaporation air temperature, (c) ice-core site condensation air temperature, and (d) ice-core site surface air temperature. All records are resampled to even 200-year resolution for visual clarity. Light shading in panels (b-d) is the absolute temperature uncertainty, while dark shading shows the relative temperature uncertainty, though the latter can be difficult to see at the scale plotted here given the larger magnitude of point-to-point variability. For more detail see Figs. A17, A18, and A19.

Discussion
Using the self-consistent, nonlinear temperature reconstruction technique for eight different ice-core sites, we next investigate the patterns of Southern Hemisphere temperature change through time. In Fig. 9 we show the nonlinear reconstructions of Antarctic surface temperature and moisture source evaporation temperature for the eight icecore records. At the WDC site in West Antarctica there is an independent estimate for the magnitude of glacialinterglacial temperature change from the borehole temperature profile . Our results are in good agreement with those findings in both the absolute value of reconstructed temperatures and the magnitude of glacialinterglacial 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,  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. Figure 10. Differences between reconstructed (a) T s and (b) T 0 using different reconstruction techniques for multiple core sites. Blue lines show the difference between our full nonlinear reconstruction and a linear reconstruction using β and γ slopes linearized around Holocene conditions. Red lines show the difference between our full nonlinear reconstruction and a linear reconstruction using slopes linearized around glacial conditions. Purple lines show the difference between the two linear reconstructions. Buizert et al. (2021) used a similar approach as  (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 , for the South Pole). We emphasize that the firn modeling approach cannot simultaneously satisfy all observational constraints, as discussed in , 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 ice-core records (Fig. 13). We weight the records equally; we do not adjust for the spatial distribution of the cores or weight by area, latitude, or elevation.
In Fig. 13 one can see that the Antarctic-wide average surface temperature change during the last deglaciation was considerably larger than the concurrent temperature change in the mean moisture evaporation source. In fact, average deglacial change in Antarctic surface temperature was about 3 times as large as the changes in evaporation temperatures, while changes in condensation temperature were about twice as large as the evaporation temperature changes.
In Fig. 14 we show the pattern of glacial-interglacial temperature change across the Antarctic continent. The magnitude of warming since the Last Glacial Maximum is calculated as the temperature difference between the late Holocene (LH, 0-4 ka) and Last Glacial Maximum (19-23 ka) for comparison with other proxy reconstructions. There may be some uncertainty in the relative magnitudes of these changes owing to offsets in the individual timescales of each record. While the relative magnitudes of interglacial change depend on the exact time periods used in the differencing, the pattern of changes in surface temperature across the continent is robust. Disentangling the full suite of causes for these temperature changes is a topic of future work but may certainly include changes in the local energy balance, heat and moisture transport, and ice sheet topography.  Tables 1 and 2). Linear methods labeled U12 for Vostok, EDC, and EDML were calculated by a simple Rayleigh-type model (Uemura et al., 2012). Reconstructions U12a-e for Dome Fuji represent a sensitivity study from Uemura et al. (2012). Reconstructions S03 and S09 are from Stenni et al. (2004Stenni et al. ( , 2010.
We find smaller glacial-interglacial temperature change for East Antarctic sites compared to previous reconstructions. Our results show that the surface temperatures of the lower, warmer areas of West Antarctica warmed significantly more than the higher, colder East Antarctic since the LGM. For example, the average surface temperature warming between the LGM and LH for the two lowest sites in our reconstruction, WDC and Siple Dome, is 11.6 • C. The average warming at the two highest sites of Dome Fuji and Vostok, however, is significantly less and just 6.9 • C or 59 % of that at the lower sites.
In Fig. 14 we plot the magnitude of warming since the LGM of Antarctic moisture source evaporation air temperatures for all ice-core records as a function of the mean latitude of the moisture source distribution for each site (based on a water-tagged general circulation model (GCM) simulations; see the Appendix). Additionally, we calculate the change in sea surface temperature (SST) during evaporation (red circles Fig. 14a) using the T 0 -SST relationship from Figure 12. Differences between our full nonlinear reconstruction and multiple previously published linear reconstructions (Stenni et al., 2004Uemura et al., 2012) of (a) ice-core site surface temperature, (b) site condensation temperature, and (c) evaporation source temperature for multiple core sites. Full reconstructions are shown in Fig. 11. our model for comparison to other SST proxy reconstructions. While plotted as points, note that these changes in moisture source temperature reflect the integrated warming over the moisture source area. Modern moisture source distributions for each site are indicated by relative histograms along the latitude axis in Fig. 14a (see Fig. 1 for further information), which are colored to correspond to the ice-core site. The moisture source points are plotted at the longitude of the respective ice-core sites in Fig. 14b, though in reality the moisture source distributions (MSDs) have a significant meridional extent (see extended data in Fig. 8, Buizert et al., 2018) that is often asymmetrically to the west owing to the westerly winds. Changes in 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 ice-core T 0 reconstructions have low spatial resolution owing to broad moisture source distributions, they benefit from the temporal resolution and precision of the ice-core age scales compared to other proxy records of temperature from the Southern Hemisphere midlatitudes.
It is clear from Fig. 14 that the Antarctic continent warmed 2 to 3 times as much as the Southern Hemisphere midlatitude moisture source areas since the Last Glacial Maximum. This result is in line with other paleoclimate reconstructions, as well as modeling of the pattern of polar amplification since the LGM Otto-Bliesner et al., 2006). In particular, our estimates of moisture source region changes agree with completely independent estimates from the MARGO compilation of SST changes (MARGO Project 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 multi-model mean pattern of Southern Hemisphere polar amplification from PIMP3 simulations is shown in Fig. 14. There is broad similarity to our reconstructions, though there are important differences as well. The spread in temperature change about the zonal mean over both the Antarctic and ocean surface is similar between the model and the reconstructions. Our re-constructions show more warming in the ice-core moisture source areas equatorward of the polar front than the PMIP3 mean and less warming over the surface of West Antarctica. We note that the magnitude and pattern of modeled Antarctic surface warming is predominately a function of imposed changes in ice sheet surface elevation to PMIP3 experiments. The extreme warming seen in parts of West Antarctica in the PMIP3 model mean (e.g., > 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 GCM-modeled 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.

Conclusions
Ice-core records of the stable isotopes of water provide detailed histories of Earth's climate. Both qualitative and quantitative interpretation of these records requires understanding Figure 14. Spatial pattern of temperature change since the Last Glacial Maximum. (a) Warming between 19-23 and 0-4 ka for ice-core site surface temperatures (colored circles with black outline corresponding to different ice-core sites as shown in map inset, uncertainty in the temperature change shown as error bars), moisture source evaporation air temperature (colored circles, red outline and uncertainty), and moisture source sea surface temperature (red circles). Moisture source warming is plotted at the mean latitude of the modern moisture source distributions for each ice-core site, while the latitudinal extent of each moisture source is indicated by the relative histograms along the x axis. Sea surface temperature warming from the MARGO compilation of SST estimates from ocean sediment cores is shown as open black circles. (b) Ice-core site surface temperature changes and moisture source sea surface temperature changes are shown as large colored circles with black outline. MARGO compilation SST changes are shown as small colored circles. (c, d) Spatial pattern of temperature change from the multi-model mean PMIP3 simulations of the LGM and pre-industrial. (c) The multi-model mean for all grid points is shown as grey dots with the zonal mean in black. Estimates of T s and T 0 from the ice-core reconstructions are shown as colored circles as in (a) for reference. 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 water-isotope variability and to disentangle the combined influences of the source and site temperatures. To date, most water-isotope temperature inversions have assumed linear relationships (Kavanaugh and 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 water-isotope models that underlie most temperature reconstructions, there are inherent nonlinearities in the isotope-temperature relationships. Ignoring these nonlinearities distorts reconstructed temperature variability. In the case of evaporation source temperature changes, these distortions may be a significant fraction of the total reconstructed variability.
There is a long-standing debate regarding the interpretation of "spatial" and "temporal" slopes in the water-isotopetemperature relationship (e.g., Jouzel et al., 1997). These discussions are conceptually useful. However, while space and time are obvious coordinates through which to understand climate, they are not the most relevant for water-isotope fractionation. Neither space nor time can independently cause water to change phase and fractionate.
The fundamental dimension through which to understand water-isotope fractionation is temperature. In this study we use a relatively simple model to investigate the relationships of water isotopes in precipitation to temperature. While the distinction between temporal and spatial slopes is not directly addressed in this context, we are able to resolve the core question: is the water-isotope-temperature relationship fixed? It is not. But these slopes are fundamentally functions of temperature, not space or time.
Our nonlinear reconstruction technique allows for the estimation of absolute temperatures in the past, in addition to their variability, and is corroborated by independent temperature constraints. By taking into account the inherent nonlinearities of water-isotope fractionation we are better able to constrain evaporation source region changes. Our reconstructions reveal a spatial pattern of temperature change across the Antarctic continent in which West Antarctica warmed significantly more than East Antarctica since the Last Glacial Maximum. Further, our reconstructions provide insight into the spatial pattern of polar amplification, suggesting that the warming since the LGM in Antarctica was 2 to 3 times that in the midlatitudes.

Appendix A: Simple Water Isotope Model
The Simple Water Isotope Model (SWIM) is based on existing numerical Rayleigh-type distillation models (Merlivat and 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 pseudo-adiabatic 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 ERA-Interim reanalysis (Dee et al., 2011). Correlations with surface air temperature are better defined than spatial correlations and give greater flexibility to the model for use in different climate states. Surface air temperature and SST are extremely well correlated in the reanalysis (Fig. A1) with a well-defined nearly linear relationship over most of the temperature range, except where the SSTs asymptote to the freezing point of seawater. The relationships are nearly identical between the NCEP and ERA reanalysis products.
Relative humidity gradients in the modern climate are fairly weak, though surface RH over the ocean is consistently higher at lower surface temperatures on climatological timescales. While variable on short timescales, RH appears to be largely invariant on timescales longer than interannual (Dai, 2006;Vimeux et al., 2002). We find that the over-ocean surface relative humidity is systematically about 5 % higher in the NCEP reanalysis compared to the ERA reanalysis, though the relationship to surface temperature is similar.
Given a specified initial air temperature, 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 • ≤ T a ≤ 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 high-order polynomial fits. All methods show effectively indistinguishable relationships in both reanalysis products. We calculate the uncertainty in these fits and test the model's sensitivity. Our base model uses the cubic spline method, which is least susceptible to edge effects and phase shifting, and it maintains a smooth first derivative.
We find some differences in the T a -to-SST 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 -to-RH 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 water-isotope 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: Figure A1. Climatological correlations between T a , SST, and RH. (a) Annual mean climatological surface air temperature, T a , and sea surface temperature, SST, from the NCEP/NCAR reanalysis (Kalnay et al., 1996). Light blue dots are Northern Hemisphere (NH) grid points, while dark blue dots are Southern Hemisphere (SH) grid points. The polynomial fit for the SH is in red. The error estimate for the fit, σ err in orange, is the standard deviation of the misfit in the model. The 1 : 1 line is shown in black. (b) Same as on the left but for the surface air temperature over ocean, T a , and surface relative humidity (RH) over oceans.
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 pseudo-adiabatically, defining a pressure trajectory, P as a function of temperature, T (Fig. 2). As the air parcel cools, moisture above saturation is removed and the latent heat released during the phase change keeps the air parcel warmer than in an otherwise equivalent isobaric pathway. Following the pseudo-adiabatic assumption, we consider no other heat sources to the air parcel, and moisture is removed immediately after condensation. Below we investigate in the influence of air parcel mixing on our model framework, which is a relaxation of the adiabatic assumption. We calculate a pseudo-adiabat following the iterative routine described in Bakhshaii and Stull (2013) but taking into account the saturated vapor pressures of both ice and liquid water condensate. The temperature-dependent saturated vapor pressures of ice and liquid water (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 R d R 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 isotope-enabled GCMs, approximate the temperature dependence of cloud ice-liquid fraction as piecewise linear functions (Hu et al., 2010), while others use smoothly varying error integrals (Ciais and 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 temperature-dependent and calculated for the liquid and ice phases individually. Effective values for the parcel as a whole are calculated from the mixing fractions of each phase (Kavanaugh and Cuffey, 2003). For example, r s(eff) = r s(ice) F (ice) + r s(liq) F (liq) , 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 (temperature-dependent) fractions of each phase of condensate.
Moisture is removed along the temperature pathway owing to the temperature-dependent changes in the saturated mix- dT . 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 water-isotope fractionation scheme used in SWIM. We model equilibrium and the kinetic fractionation of the 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(Majoube, , 1971Merlivat 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 17 O xs = δ 17 O − 0.528 × δ 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 steady-state disequilibrium, significant kinetic fractionation occurs during evaporation from the ocean. Kinetic fractionation depends both on the relative humidity and the wind speed at the air-ocean interface during evaporation (Merlivat and 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 H 18 2 O and H 16 2 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. (2008Uemura et al. ( , 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 18 α diff = 1.007 ± 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., 2008Uemura et al., , 2010Liu 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 temperature-dependent φ diff and a constant φ diff = 0.88 are small (< 1 ‰ for initial δD of vapor).
For the fractionation of 17 O 16 O we use the following relationships, which are backed by both theory and empirical observation Luz, 2005, 2007): 17 α eq = 18 α 0.529 eq for vapor and liquid in equilibrium and 17 α diff = 18 α 0.518 diff 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., , 2010Liu 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 sur-face, 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 = −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 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 v = (1−θ )R e +θ R 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  Positive values on the color scale reflect the fraction of moisture from a 5 • C warmer moisture source mixed with the local moisture source (50 = 50 % moisture from local source + 50 % moisture from 5 • C warmer-than-local moisture source), while negative values reflect the fraction of moisture from a 5 • C colder moisture source mixed with the local moisture source (−50 = 50 % moisture from local source +50 % moisture from 5 • C colderthan-local moisture source). Model results use NCEP/NCAR reanalysis for SST and RH climatology. Unmixed model results for a local closure assumption are shown as a black dashed line and a global closure assumption as a solid black line. and RH. These closure assumptions represent the bounds of a well-mixed and unmixed atmosphere or something in between. The amount of mixing in the boundary layer could change with location and with climate mean state. Rather than tying our model to one closure assumption, we view mixing at evaporation as an inherent uncertainty.
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 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 temperature-dependent α diff to reduce this misfit, we prefer a fixed α diff to avoid over fitting a relatively small dataset in the absence of further evidence or physical reasoning. While we do not consider temperature dependence of α diff , we do consider a range of α diff as an inherent uncertainty in our model and account for this in the uncertainty analysis of our temperature reconstructions as discussed in Sect. A9.
We note that it is of course also possible that other factors, rather than temperature dependence of α diff , could account for the apparent misfit, such as a difference between the ship-measured RH and 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 ship-based measurements. Given the magnitude of the misfit, it is also possible that spatial or tem-poral 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 pseudo-adiabatic cooling from the source (Dansgaard, 1964) . Thus, In general, condensation occurs in the model at saturation, and thus the temperature-dependent equilibrium fractionation factor α eq is used in Eq. (A9). However, at cold conditions there may be supersaturation of vapor over ice, leading  (Schmidt et al., 1999). (j-l) Sensitivity of modeled isotope values of vapor to a range of 18 α diff values from 1.007 to 1.010. Also shown are results using the first-principles-based diffusivities of Hellmann and Harvey (2020) (HH20, black, using n = 0.27 in Eq. A4). Local closure is assumed in all panels.
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 D 16 D 18 = 1.0285 during moisture transport, representative of pure molecular diffusion and ignoring the negligible ventilation effect. Likewise we use D 1 D 2 = 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 temperature-dependent φ diff used in the evaporation scheme. We prefer this value for simplicity, consistency with earlier work, and lack of experimental measurements of φ diff at the colder temperatures experienced during transport (Luz et al., 2009). However, to investigate the effect of assuming a constant φ diff , we also use a configuration of the model with the temperature-dependent diffusivity ratios of Hellmann and Harvey (2020) based on molecular kinetic theory. We can achieve essentially identical results using this physically based temperature dependence of the diffusivity ratios and assuming no temperature dependence, with exceedingly minor changes to the supersaturation function (see Sects. A2.3 and A4); a change of b in the supersaturation function of less than 0.00015 is sufficient to negate any difference between these assumptions. We thus consider the effect of temperature dependence in these diffusivities to be subsumed within the uncertainty associated with tuning the supersaturation function (see Sects. A4 and A9) In the mixed-phase portion of the transport pathway, the effective fractionation factors are determined by the mixing fractions of ice and liquid condensate. Following Kavanaugh and Cuffey (2003), The temperature dependence of the ice and liquid fraction is shown in Fig. A7a and based on satellite observations (Hu et al., 2010).

A2.3 Supersaturation
The supersaturation of vapor over ice is a critical parameterization in water-isotope distillation models. The true relationship of supersaturation to environmental conditions is the result of complex cloud microphysics (Hong et al., 2004). Because of its strong influence on water-isotope fractionation and the uncertainty in the underlying physics, the supersaturation is often parameterized to depend on temperature and tuned to fit water-isotope models to observations (Jouzel and 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 water-isotope data could not distinguish among possible functional forms of the parameterization (e.g., linear, exponential). Their linear parameterization has been used extensively in water-isotope models since: where a and b are tuned to fit observational data. It is important to note here that the prescribed mixing of liquid and ice in the cloud implies a supersaturation of vapor over ice that follows the blue curve shown in Fig. A7b, which is inconsistent with the supersaturation driving kinetic fractionation as prescribed in Eq. (A14). The presence of both liquid and ice phases in a cloud is not the only potential source of supersaturation. The lack of condensation nuclei, for example, allows supersaturation to remain high in cold, ice-only conditions (Hong et al., 2004), rather than returning to unity as the cloud becomes entirely ice-phase. It is common for water-isotope models, even those in some GCMs (e.g., Schoenemann et al., 2014), to have multiple variables equivalent to supersaturation in different aspects of the same model, such as the isotope fractionation and precipitation schemes, which may not be self-consistent. Because the environmental supersaturation experienced by the air parcel is related to the relationship between temperature and moisture removal (that is d ln(f ) dT ) and the supersaturation driving kinetic fractionation relates temperature to δ values (largely through Eq. A12), an inconsistency in the model's view of supersaturation can influence the modeled 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-water-isotope relationship that results if inconsistent relationships between saturation and temperature are used in the precipitation and fractionation schemes, which is generally incompatible with observations. In line with previous work (Jouzel and 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 non-supersaturated conditions in the ice-only portion of the cooling pathway, the simultaneous transition to equilibrium-only fractionation would drive a slope of ∂δ 18 O ∂δ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 ice-core sites

A3.1 Moisture source distributions
The atmospheric circulation transports moisture poleward of ≈ 30 • S (Fig. 1a). The mean evaporation latitude of moisture that precipitates at any given site can be estimated from moisture-tagged GCM experiments (Markle et al., 2017). The difference between the mean latitudes of moisture evaporation and precipitation steadily increases between the sub- × T . Note that the atmosphere is subsaturated with respect to ice at temperatures above 0 • C but that our pathways do not include ice at these temperatures. tropics and the pole ( Fig. 1a and b). The mean evaporative latitude of moisture that precipitates in Antarctica ranges between 44 and 50 • S (Fig. 1c). The surface elevation of the ice sheet is a strong predictor of the mean latitude of precipitation, with higher-elevation sites having more equatorward moisture sources (Fig. 1d) due to topographic isolation (Sodemann and Stohl, 2009;Bailey et al., 2019). There are some notable asymmetries in this general pattern. The large embayments are areas of comparatively high-latitude mean moisture source, such as the Victoria Land coast in the Ross Sea region.
The mean moisture source latitude is, however, not the full story. The moisture reaching any Antarctic site does not originate from just a single mean source latitude or follow a single temperature pathway. The contribution of moisture evaporated from different latitudes to the final precipitation at a site defines a moisture source distribution (Markle et al., 2017), which reflects the combination of the spatial pattern of evaporation, cumulative rainout, and the influence of atmospheric circulation. Here we diagnose annual mean moisture source distributions (MSDs) as a function of latitude from a moisture-tagged run of the Community Atmosphere Model (CAM) for East and West Antarctic sites (details are given in Markle et al., 2017), as shown in Fig. 1e. Moisture source distributions derived from other methods like trajectory modeling are similar (e.g., Sodemann and 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 ice-core sites is less than 10 • of latitude, while the moisture source distribution for any one site spans over 40 • of latitude. Local evaporation is a small contribution to the moisture precipitating at Antarctic sites. On average moisture is transported more than 20 • of latitude to reach Antarctica.
Given the broad range of evaporative conditions that contribute to moisture precipitating at an ice-core site, what is the meaning of the T source that can be reconstructed from water-isotope records? It is the moisture-weighted evaporative temperature, determined by the convolution of the spatial pattern of the MSD and the underlying surface temperatures (Fig. A28, Markle et al., 2017). Both surface temperatures at fixed locations and the pattern of the MSD can change independently. The water-isotope records alone do not allow the disentanglement of these two patterns, which may have different temporal evolution (Markle et al., 2017).
To understand the moisture transport and water-isotope distillation to Antarctic sites it is important to consider evaporation from the range of conditions comprising the moisture source distribution. We thus use an ensemble of temperature pathways for Antarctic precipitation defined by a range of Antarctic condensation temperatures as well as the broad range of evaporation temperatures underlying the Antarctic moisture source distributions. The means of these distributions vary across the continent.

A3.2 Condensation site conditions
During transport, moisture is cooled from initial surface air temperature at evaporation to subsequent condensation temperatures. The condensation temperature is not the same as the surface temperature where that precipitation falls. Indeed, there is a vertical and temporal distribution of condensation contributing to precipitation that falls at any point on the surface, analogous to the horizontal and temporal distribution of evaporation contributing the moisture ultimately transported to any precipitation site. What is the meaning of T c reconstructed from ice-core records? It is the vertical profile of temperature weighted by the vertical profile of condensation that yields net accumulation to a site. The weighted condensation temperature has a distinct relationship to the surface temperature across the globe.
Antarctica has strong climatological inversions such that temperature aloft is often warmer than the surface (Connolley, 1996). Masson-Delmotte et al. (2008) review the rela-tionship 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 ERA-40 reanalysis  and MAR, a high-resolution mesoscale model forced by ERA-40 (see Fig. 8, Masson-Delmotte et al., 2008). In both models the upper bound of the Antarctic condensation temperature appears to be set by the peak inversion temperature, though condensation temperatures are on average colder than the peak inversion temperature (meaning simply that condensation occurs at a range of temperature up to the peak inversion temperature). Masson-Delmotte et al. (2008) calculate a best fit of the surface to condensation temperature slope of 0.65 • C • C −1 in the ERA-40 data, consistent with previous work that found a slope of 0.67 • C • C −1 (Connolley, 1996;Jouzel and Merlivat, 1984). The spread of condensation temperatures in the higher-resolution MAR model suggests colder condensation temperatures than in the lowerresolution reanalysis . The strength of the Antarctic inversion diminishes with increasing surface temperature (Connolley, 1996), and relatively warm Antarctic surface temperatures (e.g., > −20 • C) are associated with condensation temperatures colder than the surface temperature .
We analyze the relationship between surface temperature and condensation temperature in monthly MERRA-2 reanalysis from 2008 through 2017 (Gelaro et al., 2017). We show the climatological zonal-mean vertical profile of air temperature, the sum of the convective and large-scale precipitation source production rate, and the rate of reevaporation and sublimation of precipitation in Fig. A8. The relationship between the climatological weighted condensation temperature and the surface air temperature at every grid point is shown in Fig. A9. Note that this calculation accounts for the seasonality of precipitation throughout the atmospheric column, as well as the reevaporation and sublimation of falling precipitation. Ignoring reevaporation and sublimation leads to qualitatively similar results.
The primary take-away is that the MERRA-2 data show a generally linear relationship between condensation and surface temperature for typical Antarctic surface temperatures. That relationship, however, is not linear at warmer surface temperatures. Indeed, even at surface temperatures below zero, the relationship is not strictly linear, but rather steepens with decreasing temperature. The relationship between the surface air temperature and the weighted condensation temperature (for surface temperatures below −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-0.75 • C • 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 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 MERRA-2 dataset, even to some extent at warm surface temperatures. This is despite the potential for offsets between the reported surface temperature and that at the time of the precipitation for the MD08 and GNIP precipitation samples, differing time intervals, potential moisture biases in the column in MERRA-2 (Bosilovich et al., 2017), and the lack of processes in our isotope distillation model that should be important, for example, in tropical convection or that might alter Antarctic precipitation after deposition.
Examining all modeled condensation temperatures for samples in the MD08 and GNIP datasets with reported sur- face temperatures below 15 • C, we find slopes between 0.62 and 0.67 • C • C −1 . For just the Antarctic precipitation samples in the MD08 dataset we find a best-fit slope between the reported surface temperature and our modeled condensation temperature of 0.67-0.69 • C • C −1 (Fig. A9), depending on the model assumptions (0.69 • C • C −1 under our base assumptions) and whether below-freezing source evaporation is included (see below), in good agreement with previous Antarctic observational studies (Connolley, 1996;Jouzel and Merlivat, 1984). These slopes sit well within the range found in the MERRA-2 data. We also reconstruct the condensation temperatures for the topmost samples from several deep ice cores and compare those to the reported annual average temperatures for those sites (Fig. A9). We find a best-fit slope between 0.68 and 0.70 • C • C −1 , depending on whether we average samples from the last 50 or 100 years, though only five points describe these lines.
Based on the above results we use the equation T c = 0.69 • C • C T s −8.2 • C as our base estimate to reconstruct Antarctic surface temperatures,; however, we consider an uncertainty of ±0.02 • C • 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 temperature-dependent isotope slopes in Fig. A11.

A3.3 Seasonality
Does seasonality in the hydrological cycle systematically bias climatological information in ice-core 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 water-isotope ratios of precipitation. As discussed above, the critical comparison is between the annual mean surface temperature and the condensation-weighted temperature, integrated over both time and altitude. Our analysis of the MERRA-2 reanalysis data (Fig. A9) takes seasonal variation in precipitation and the vertical temperature profile into account. Differences in seasonality of condensation at different sites contributes to the spread around the central relationship. It is nevertheless useful to investigate potential bias in the 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 ERA-Interim reanalysis as well as global precipitation in the MERRA-2 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 month-to-month and year-to-year variance in snowfall is large compared to the climatological seasonal cycle. The stochastic sampling of the seasonal cycle in surface temperature overwhelms the potential bias introduced by the average seasonality of precipitation. Further, the timing of the climatological annual cycle in snowfall varies across the continent, whereas the annual temperature cycle is quite coherent. The potential for seasonal bias thus varies dramat- ically between sites, even in the absence of the dominant stochastic sampling.
We compare the annual average surface temperature to the precipitation-weighted annual temperature at every grid point in the Southern Hemisphere. The mean bias for sites with typical Antarctic surface temperatures is less than 0.33 • C, with the precipitation-weighted temperature being slightly colder on average. We find no systematic dependence of this bias on the surface temperature itself. While individual sites do show differences up to 4 • C over the interval examined, our analysis does not suggest that such differences are persistent at a site. None of these analyses of monthly data take into account potential biases at the scale of individual precipitation events. The intermittency of Antarctic snowfall likely complicates the relationship between condensation-weighted and annual mean temperature at the seasonal and annual scale. At the same time, however, precipitation intermittency reduces potential seasonal biasing at climatological timescales by degrading any coherence in the seasonal cycles of accumulation and temperature.
Could seasonality in evaporation bias reconstructed source region T 0 ? We examine the seasonality of Southern Hemisphere evaporation in the monthly MERRA-2 reanalysis by comparing the annual average over-sea surface temperatures to the evaporation-weighted annual temperatures. Between 35 and 65 • S, the bulk of Antarctic moisture sources and evaporation-weighted surface temperatures are slightly colder than mean annual surface temperatures (the mean difference is 0.123 • C, with 95 % of points between −0.25 and 0.5 • C). South of the climatological sea ice zone, mean evaporation temperatures are a couple degrees warmer than mean annual surface temperatures on average, though our moisture tagging analysis suggests that these areas contribute at most a couple percent of the total moisture arriving at typical Antarctic sites.

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 i = a + b × T , set a = 1, and tune the slope, b. The supersaturation has a strong influence on the kinetic fractionation (Eq. A12) and thus the relationship between δD and δ 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 water-isotope measurements of precipitation and surface snow from Antarctica and around the globe. The bulk of this compilation is that published by Masson-Delmotte et al. (2008). We include additional published surface snow and precipitation measurements from the GNIP database (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 L2120-i analyzer). 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 δ 18 O in a global dataset for precipitation. They use a second-order polynomial fit, which is the basis for the logarithmic deuterium excess parameter (Uemura et al., 2012;Markle et al., 2017) (Eq. 4). From Eq. (A9), we can see that the theoretical relationship between δD and δ 18 O, given any amount of distillation, depends on the ratio of exp( D α tot ) exp( 18 α tot ) , where each α tot is itself a nonlinear function of temperature as outlined throughout Sect. A2.2. The ratio of exponential functions can be estimated to any arbitrary degree of accuracy with a polynomial function.
Our modern dataset includes several new sets of measurements in addition to those used in Uemura et al. (2012). We find similar coefficients in a second-order polynomial fit between δ 18 O and δ D in our larger dataset compared to those found by Uemura et al. (2012): A = −29.2 and B = 8.45 (see Eq. 4). Because these coefficients are not significantly different than those previously published and for consistency with that work, we use the coefficients found by Uemura et al. (2012) (A = −28.5 and B = 8.47) to define a logarithmic deuterium excess parameter, d ln . We find no benefit or justification for using higher-order fits to this dataset.
We run SWIM to produce an ensemble of temperature trajectories representing a wide range of possible evaporation and condensation temperatures (T 0 varies from 0 to 28 • C; T c from 27 to −60 • C). We then compare the resulting cloud of modeled δ 18 O, δ D, and d ln , finding a second-order polynomial fit between the modeled δ D and δ 18 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 δ 18 O (Eq. 4). This is easily visualized in a plot of δ 18 O and d ln (e.g., Fig. A12a), as the average δ 18 O-to-d ln relationship is flat in measured samples by definition.
Using the local closure assumption we find an optimal tuning of S i = 1 − b × 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 δ 18 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 second-order polynomial fit between modeled δ D and δ 18 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 higher-elevation, colder Antarctic sites likely have more equatorward MSDs (Fig. 1), we should expect more depleted δ 18 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 δ 18 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 δ 18 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., 0.005 ≤ b ≤ 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, temperature-dependent parameterization of supersaturation is both simple and widely used, the physical processes determining supersaturation are complex. To understand the sensitivity of the model to this parameterization we also test a nonlinear parameterization of supersatu- 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 second-order term, c, those of order 5 × 10 −6 • C −2 , are reconcilable with the observed δ 18 O-to-d 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 moisture-weighted average of a set of independent pseudo-adiabatic pathways, defined by a range of evaporation and condensation temperatures, can approximate the conditions experienced by precipitation falling at a polar site. We now aim to test the limits of this approximation and assess the influence of atmospheric mixing on the isotopic composition of air masses within the simple model framework.
The influence of air mass mixing during evaporation is considered in the discussion of the closure assumption above. Here we consider mixing during transport of air masses with different initial evaporation conditions, different condensation histories, and different temperature, moisture content, and isotopic values at the time of mixing. The central question is whether the processes of mixing can result in isotopic values of final precipitation that are significantly different than the moisture-weighted average of precipitation from unmixed pathways. There are three processes associated with mixing to consider.
1. Non-uniqueness. Parcels that experienced different evaporation conditions can arrive at a condensation site of a given temperature with different isotopic values.

2.
Mixing-induced condensation. Mixing two saturated or undersaturated air parcels of different temperatures may result in an oversaturated mixed parcel due to nonlinearity in the Clausius-Clapeyron relationship. This process leads to additional condensation and fractionation (as well as warming due to latent heat release) and thus a more depleted isotopic value for a given temperature compared to the moisture-weighted average of unmixed pathways.
3. Nonlinear mixing. For isobaric mixing of equal-massed air parcels, the final mixed temperature reflects their mass-weighted average (plus the effect of any latent heat release). However, the relative abundances of water isotopes mix with the moisture content of the air parcels rather than their total mass. Thus, while the temperatures have mixed linearly, the isotopic values of the resultant mixed parcel will be weighted nonlinearly toward those of the warmer and wetter parcel. The result is that the mixed parcel has a less depleted isotopic value for a given final temperature compared to the moistureweighted average of the unmixed parcels.
Moisture-weighted differences between mixed and unmixed parcels only occur when air masses of different temperatures mix. Physically this may represent colliding fronts at synoptic scales. We consider two air parcels evaporated from identical starting conditions but which mix at different temperatures. The influence of processes 2 and 3 is largest when the relative humidities of both parcels are at saturation, and the magnitude of the influences increases both as the difference in temperature between the two parcels increases and as the absolute temperature of the parcels increases. As temperature decreases the nonlinearity of the mixing of moisture approaches the linear (mass-weighted) mixing of temperature.
To assess the range of temperature differences associated with synoptic scales in the Southern Hemisphere, we examine the difference in daily mean 2 m air temperature from the ERA-Interim reanalysis (Dee et al., 2011) over the Southern Hemisphere oceans. In summer, the day-to-day temperature differences have a standard deviation of less than 0.9 • C, and in winter they have a standard deviation of less than 1.5 • C (other reasonable metrics of synoptic-scale variability such as lagged 2 or 5 d temperature differences or grid-point-to grid-point differences are similar). While there is surely the potential for strong mixing for any given synoptic event, for the purposes of paleoclimate reconstruction we are interested in the long-term average of many storm events and thus the statistics of mixing generally.
We use the simple model to assess the range of final isotopic values of precipitation that can arise from two air parcels, which evaporated at the same initial conditions, then mixed at a range of different temperatures during transport. Air parcels are evaporated at specified initial air temperature (T 0 = 10 • C), cooled, and randomly mixed at any com-bination 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 moisture-weighted means of the mixed distributions are shown as vertical black lines, while the unmixed final values are shown as vertical red lines. The means of the distribution are less depleted than the unmixed parcel owing largely to process 3. The influence of process 2 can also be seen in the additional peak at more depleted values. These distributions vary as a function of both T 0 and T c , though the differences in moisture-weighted means between mixed and unmixed parcels are relatively small and consistent.
The idealized tests above show the influence of mixing (at different temperatures) of air parcels that evaporate from identical source conditions. Perhaps more realistically, air parcels from different sources can mix at different air temperatures, in which case all three mixing processes above are important. This can act to broaden the distribution associated with a given moisture-weighted mean isotopic value. In Fig. A14e-h we show the distribution of the isotopic value of precipitation at −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 moisture-weighted 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 moisture-weighted means of mixed and unmixed parcels are actually smaller. This is presumably because the skewed influence of process 3 (which drives the persistent bias above) contributes less to the total distribution. The relationship between δ 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 ice-core 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 200-year resolution for visual clarity. The light shading around the reconstructions represents the uncertainty in absolute temperature, while the darker shading represents the uncertainty in relative temperature variability as described in Sect. A9. Note that in some cases sample-to-sample variability is larger than the relative uncertainty, making the shading difficult to see at the scale plotted here.

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). We also show the difference in the reconstructions using the Holocene and glacial linearizations.
We compare the nonlinear temperature reconstructions of all eight ice-core sites to linearized reconstructions using SWIM results for Holocene conditions as described in Sect. 4.3. 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). Re-constructed 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 ∂d xs ∂T 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  Figure A17. Source region evaporation air temperature reconstruction for all ice-core sites. The light shading around the reconstructions represents the uncertainty in absolute temperature, while the darker shading represents the uncertainty in relative temperature variability as described in Sect. A9. Figure A18. Condensation air temperature reconstruction for all ice-core sites. The light shading around the reconstructions represents the uncertainty in absolute temperature, while the darker shading represents the uncertainty in relative temperature variability as described in Sect. A9.
ice-core 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 Figure A19. Surface air temperature reconstruction for all ice-core sites. The light shading around the reconstructions represents the uncertainty in absolute temperature, while the darker shading represents the uncertainty in relative temperature variability as described in Sect. A9. 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 α 18 diff during evaporation is important in setting the initial isotopic values of vapor. While we find the value 1.009 to give the best fit to modern observations, values within a small range may be defensible (Fig. A6). The local closure assumption used in the evaporation scheme has known limitations (Risi et al., 2010), representing an end-member scenario for possible evaporative conditions. While less applicable to past climate mean states, the global closure assumption provides an extreme test of the model's sensitivity. Using the global rather than local closure assumption can lead to differences in reconstructed absolute 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 deter- mine an initial relative humidity given an initial air temperature; colder surface air temperatures over the ocean are associated with slightly higher relative humidity. We show the difference in reconstructions due to using either the NCEP or ERA-Interim reanalysis. We can also ask how different our reconstructed 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 ran-  dom pairs of pseudo-adiabatic 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 moisture-weighted 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 tech-  nique 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 T c = T c,M − T c and T 0 = T 0,M − T 0 . We test a range of mixing and threshold values in multiple Monte Carlo simulations. In all cases the mean values of T c,M − T c and T 0,M − T 0 are very near zero (< | ± 0.06 • C| for T c , and < | ± 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 pseudo-adiabatic cooling pathways. Including moisture sources with T 0 < 0 • 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 0,M − T 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 re- Figure A25. Temperature reconstruction sensitivity to model parameterizations. T c and T 0 reconstructions using WDC δ 18 O and d ln (resampled to 10-year resolution) as well as SWIM runs. Black lines show the base model in all panels, while colored lines are results in which model parameters and assumptions are varied: (a-b) evaporation condition correlations based on NCEP (black) and ERA (red) reanalysis data; (c-d) global (blue) and local (black) closure assumption during evaporation. (e-f) A range of values for 18 α diff ; (g-h) a range of values for the b parameter in the supersaturation parameterization, as well as a nonlinear parameterization (c) as described in the text; (i-j) several versions of the precipitation parameterization in which all moisture is removed above saturation (sat), all moisture is removed above initial RH (RH = RH 0 , constant RH along path), and all moisture is removed above fixed 80 % or 90 % RH.
constructions 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 18 a diff = 1.009±0.001. We estimate the uncertainty introduced by the precipitation parameterization as the mean absolute difference from the base scenario of reconstructions using each of the alternate assumptions outlined in Sect. A1.2, applied symmetrically to the base scenario. We estimate the uncertainty arising from our assumed relationship between T 0 and RH 0 as the mean absolute difference in reconstructions using climatological fits from the NCEP/NCAR and ERA-Interim reanalysis.
Based on the tests in Sect. A2.1, a conservative estimate of the uncertainty arising from mixing at the evaporation source Figure A26. The influence of air parcel mixing on modeled isotope state spaces. (a, e) The water-isotope ratios of precipitation resulting from 7.5 × 10 5 stochastically mixed distillation pathways, colored by T c and T 0 , respectively. (b, f) The water-isotope ratios of precipitation resulting from the corresponding 1.5 × 10 4 unmixed distillation pathways, colored by T c and T 0 , respectively. (c, g) The difference in T c and T 0 , respectively, between the mixed and unmixed pathways as a function of the δ 18 O and d ln space. (d, h) Histograms of the differences, T c and T 0 , from all points in the δ 18 O and d ln space, resulting from mixing.
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 -to-T s relationship outlined in Sect. A3.2. We use the mean absolute difference of reconstructions using T c ∝ 0.69±0.02 • C • C T 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 ice-core 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 water-tagged GCM experiments (Fig. A28). In fact, 85 %-95 % of all moisture that arrives at Antarctic sites in our modeling initially evaporates from locations with an- Figure A27. Sensitivity of the T s reconstruction for WDC to the relationship between T c and T s . Our base scenario is shown in black (T c ∝ 0.69 • C • C T s ), while the spread associated with a ±0.02 • C • C uncertainty in the scaling factor is shown in red and gold. For comparison the results from a weaker slope of 0.65 • C • C are shown in blue. This range of scaling factors has little impact on the reconstructed temperature history and may be difficult to see. nual 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 sur- Figure A28. Annual mean temperature of moisture source regions from moisture-tagged GCM experiments. (a) Annual mean zonal-mean surface air temperature vs. latitude (blue) overlain with annual mean MSDs for several ice-core sites, colored by MSD-weighted annual mean air temperature. (b) MSD-weighted annual mean air temperature for every model grid point over the Antarctic. Note that this is not strictly the evaporation-weighted air temperature; it is weighted by location but not necessarily by time and may thus be biased to be too low. Figure A29. Isotope model results colored by initial evaporation air temperature T 0 ( • C) using base model assumptions. For initial evaporation air temperatures below 0 • C there are non-unique results in the δ 18 O and d ln space. face is shown in Fig. A29. Given a lack of isotopic vapor measurements for evaporation air temperatures much below 0 • C and that our evaporation scheme is not well calibrated for such conditions, these results are purely illustrative.
Those caveats aside, we investigate the sensitivity of our temperature reconstructions to this potential non-uniqueness in the water-isotope state spaces. We run SWIM through a large field of T 0 values from −28 to 28 • C. We can in principle resolve the non-uniqueness problem by combining reconstructions from either side of the folded surface based on the contribution of total moisture represented by each pair of non-unique paths. Knowing that below-zero moisture sources contribute far less to the total moisture reaching Antarctic sites (Fig. A28), we examine two reasonable methods of moisture-weighting the reconstructions. In the first we simply weight each pair of reconstructed temperatures by the final mixing ratio (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 higher-elevation, colder sites compared to our GCM-based 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 non-uniqueness 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 Figure A30. Temperature reconstructions for WDC accounting for below-zero evaporation air temperatures. (a) Reconstructions of T c under base conditions (blue; no contribution from T 0 < 0 • C), with the contribution from T 0 < 0 • C weighted by the final r s (red) and with a 90 % contribution from T 0 > 0 • C as well as 10 % from T 0 < 0 • C (gold). (b) Same as in (a) but for reconstructions of T 0 . using either method and less than 0.19 • C for T 0 . Attempting to account for this non-uniqueness does, however, lead to persistent mean offsets in absolute temperature; in particular, we find colder absolute values of reconstructed T 0 for all ice-core sites.
While these results are interesting, this attempt to account for non-uniqueness likely does not actually improve the absolute temperature reconstructions. Given the shape of the folded temperature surfaces in the modeled δ 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. Near-surface relative humidity in particular is not well constrained by our climatological correlations in these circumstances. Our model likely underestimates the depletion of the initial evaporate in these circumstances, meaning that the reconstruction solves for a very cold 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 con-servative estimate of the influence of non-uniqueness on our temperature reconstructions and their relative variability; the real effect is likely much smaller though difficult to quantify precisely.

A10 Comparison to previous reconstructions
We next reconstruct site and source temperatures for four East Antarctic ice-core records and compare to previously published linear reconstructions. We use records of δ 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, EPICA Dronning Maud Land (EDML) , and Dome Fuji records (Uemura et al., 2012). After seawater correction, we use the ice-core δ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 water-isotope models used for the estimation of scaling relationships. In many cases, the previously published linear inversions overestimate changes in both site and source temperature compared to the nonlinear reconstruction.
The overestimation of reconstructed temperature change by the linear reconstruction makes physical sense. The largest source of nonlinearities in the water-isotope to temperature relationships is the deuterium excess parameters, ∂d xs ∂T c and ∂d xs ∂T 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 relation-ships 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 Three-parameter 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, 17 O xs = δ 17 O − 0.528δ 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 three-dimensional 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 non-uniqueness to address; that is, a position in the three-parameter 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 , water-isotope diffusion lengths (Johnsen et al., 2000), or the concentration of aerosols , may be added to the model. These additional proxies can allow for the reconstruction of additional independent variables and the relaxation of assumptions. Alternatively it may be possible to use the same approach to optimize model parameters like the supersaturation. We leave this task to future work.
Code and data availability. The temperature reconstructions, underlying data, and Simple Water Isotope Model code are available through a public GitHub repository (doi: https://doi.org/10.5281/zenodo.6510097, Markle and Steig, 2022; https://github.com/bradley-markle/simple_water_isotope_model, Markle, 2022) and the United States Antarctic Research Program (USAP) Data Center. The ice-core data are all already available through links provided in the primary papers cited for each dataset.
Author contributions. BRM conceived of the study and, with EJS, designed the model and analyses. BRM and EJS wrote the paper.
Competing interests. The contact author has declared that neither they nor their co-author has any competing interests.
Disclaimer. Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.