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

1 Introduction

Stable-isotope ratios of water have been the foundational proxy for polar paleoclimate research for more than a half-century (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 (Masson-Delmotte et al., 2006; EPICA Community Members, 2006; WAIS Divide Project Members, 2013, 2015). Both oxygen and hydrogen have stable isotopes whose ratios (${}^{\mathrm{18}}\mathrm{O}/{}^{\mathrm{16}}\mathrm{O}$ and ${}^{\mathrm{2}}\mathrm{H}/{}^{\mathrm{1}}\mathrm{H}$) are commonly expressed as deviations, δ¹⁸O and δD, from Vienna Standard Mean Ocean Water (VSMOW) in per mill (‰):

$$\begin{array}{}\text{(1)}& \mathit{\delta }={\displaystyle \frac{R_x-R_{\mathrm{std}}}{R_{\mathrm{std}}}},\end{array}$$

where R_x is the ratio in the sample and R_std is the ratio in VSMOW.

Poleward transport of moisture by the climate system, the progressive removal of moisture from the atmosphere by condensation and precipitation, and the fractionation of water-isotope ratios during phase changes are all processes inherently linked to temperature and together underpin the use of water-isotope ratios in polar precipitation as a temperature proxy (Craig, 1961; Epstein et al., 1963; Dansgaard, 1964; Gonfiantini, 1965). The strong empirical correlation between the water-isotope ratios in precipitation and surface temperature supports this interpretation (Petit et al., 1999; Jouzel et al., 1997; Masson-Delmotte et al., 2008). Air temperatures during condensation (Petit et al., 1999; Jouzel et al., 1997) and during initial moisture evaporation (Vimeux et al., 2002) can be reconstructed from ice-core water-isotope records if the relevant scaling relationships can be determined from theory, models, or observations (Vimeux et al., 2002; Kavanaugh and 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.

1.1 Temperature reconstructions

Any interpretation of water-isotope ratios as a proxy for temperature requires a model, whether conceptual, statistical, or numerical. A conceptual model of progressive distillation and integrated fractionation (e.g., 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 δ¹⁸O and δD, d_xs=δD$-\mathrm{8}\times {\mathit{\delta }}^{\mathrm{18}}$O, and is commonly used to quantify these differences (Dansgaard, 1964; Merlivat and Jouzel, 1979).

Changes in water-isotope parameters measured in precipitation at an ice-core site, Δδ¹⁸O and Δd_xs, can be conceptualized as driven by changes in site and evaporation source temperature, ΔT_site and ΔT_source:

$$\begin{array}{}\text{(2)}& {\displaystyle }\mathrm{\Delta }{\mathit{\delta }}^{\mathrm{18}}\mathrm{O}={\mathit{\gamma }}_{\mathrm{1}}\mathrm{\Delta }{T}_{\mathrm{site}}+{\mathit{\gamma }}_{\mathrm{2}}\mathrm{\Delta }{T}_{\mathrm{source}},\text{(3)}& {\displaystyle }\mathrm{\Delta }{d}_{\mathrm{xs}}={\mathit{\beta }}_{\mathrm{1}}\mathrm{\Delta }{T}_{\mathrm{site}}+{\mathit{\beta }}_{\mathrm{2}}\mathrm{\Delta }{T}_{\mathrm{source}},\end{array}$$

where β and γ are the partial derivatives of δ¹⁸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 δ¹⁸O and d_xs (Vimeux et al., 2002; Stenni et al., 2010; Uemura et al., 2012).

1.2 Nonlinearities in isotope fractionation and the deuterium excess definition

The temperature reconstruction approach described above depends on the assumption that the parameters, β and γ, are fixed in time and independent of temperature. However, the β and γ parameters, as diagnosed from model simulations, are found to be different for different 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 δ¹⁸O and δD is not fundamental (Craig, 1961); equilibrium fractionation alone drives a nonlinear relationship between δ¹⁸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 $\mathit{\delta }{{}^{\prime }}^{\mathrm{18}}\mathrm{O}$ and δ^′D data, where ${\mathit{\delta }}_x{}^{\prime }=\ln \left(\mathrm{1}+{\mathit{\delta }}_x\right)$, and defined a phenomenological, nonlinear deuterium excess parameter:

$$\begin{array}{}\text{(4)}& d_{\ln }=\mathit{\delta }{}^{\prime }\mathrm{D}-\left(A\times (\mathit{\delta }{{}^{\prime }}^{\mathrm{18}}\mathrm{O}{)}^{\mathrm{2}}+B\times \mathit{\delta }{{}^{\prime }}^{\mathrm{18}}\mathrm{O}\right),\end{array}$$

with coefficients $A=-\mathrm{28.5}$ and B=8.47 (note that the coefficients and δ^′ values are unitless; for example, $\mathit{\delta }{{}^{\prime }}^{\mathrm{18}}\mathrm{O}=\mathrm{0.040}$ not 40, with the ‰ ignored).

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

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

2 A (relatively) simple water-isotope model

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

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

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

Our model distills moisture down thermodynamic pathways defined by temperature and pressure. Changes in water-isotope ratios are driven neither by changes in space nor time but by changes in the thermodynamic variables that cause the water to change phase. The temperature gradient of the pathway is prescribed from an initial evaporation temperature, T₀, 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₀, to other initial conditions including sea surface temperature, SST₀, and relative humidity, RH₀ (see Sect. A1.1).

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

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

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

3 Temperature-dependent slopes

To investigate the sensitivity of water-isotope ratios of Antarctic precipitation to site and source conditions, we use SWIM to model isotopic state spaces. We run SWIM through a large ensemble of temperature pathways defined by $T=T_{\mathrm{0}},T_{\mathrm{0}}-\mathrm{d}T,\mathrm{\dots },T_{\mathrm{c}}$, with dT=0.1 ^∘C. We run nearly 24 000 trajectories to fill the plausible parameter space of 0 ^∘C $\le T_{\mathrm{0}}\le \mathrm{28}$ ^∘C and −70 ^∘C $\le T_{\mathrm{c}}\le \mathrm{10}$ ^∘C. We first examine the model with base assumptions and parameterizations². 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₀, respectively, and whose z (color) dimension is the isotopic value of final precipitation in Antarctica: δ¹⁸O in Fig. 3a, δD in Fig. 3b, d_xs in Fig. 3c, and d_ln in Fig. 3d.

Figure 3 Water-isotope state spaces as a function of the boundary conditions T_c, the condensation temperature, and T₀, the mean evaporation temperature. Surface shading and contour lines are the water-isotope values of precipitation (in units of ‰) (a) δ¹⁸O, (b) δD, (c) linear deuterium excess d_xs, and (d) nonlinear deuterium excess d_ln.

The gradients of both the δ¹⁸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 δ¹⁸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₀, due to its influence on the total distillation gradient). Further, the partial slopes of δ¹⁸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₀ and T_c are shown explicitly in Figs. 4 and 5. It is important to recognize that the partial derivatives with respect to T₀ are not for the initial vapor at the point of evaporation, but for the precipitation after the vapor has passed through the distillation pathway to the final precipitation site. The sensitivity of isotopic values of precipitation to source region conditions is a function of the total distillation that the moisture experiences.

Figure 4 Partial derivatives of isotope state spaces with respect to condensation temperature: (a) $\frac{\partial {\mathit{\delta }}^{\mathrm{18}}\mathrm{O}}{\partial T_{\mathrm{c}}}$, (b) $\frac{\partial \mathit{\delta }\mathrm{D}}{\partial T_{\mathrm{c}}}$, (c) $\frac{\partial d_{\mathrm{xs}}}{\partial T_{\mathrm{c}}}$, (d) $\frac{\partial d_{\ln }}{\partial T_{\mathrm{c}}}$. Shading and contours in all panels represent the slope in ‰ ^∘C⁻¹.

Figure 5 Partial derivatives of isotope state spaces with respect to evaporation temperature: (a) $\frac{\partial {\mathit{\delta }}^{\mathrm{18}}\mathrm{O}}{\partial T_{\mathrm{0}}}$, (b) $\frac{\partial \mathit{\delta }\mathrm{D}}{\partial T_{\mathrm{0}}}$, (c) $\frac{\partial d_{\mathrm{xs}}}{\partial T_{\mathrm{0}}}$, (d) $\frac{\partial d_{\ln }}{\partial T_{\mathrm{0}}}$. Shading and contours in all panels represent the slope in ‰ ^∘C⁻¹.

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 δ¹⁸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 γ₁ and γ₂ in Eq. (2) are comparable to the slopes $\frac{\partial {\mathit{\delta }}^{\mathrm{18}}\mathrm{O}}{\partial T_{\mathrm{c}}}$ and $\frac{\partial {\mathit{\delta }}^{\mathrm{18}}\mathrm{O}}{\partial T_{\mathrm{0}}}$ in Figs. 4 and 5. Although the slope of δ¹⁸O along the condensation temperature axis, $\frac{\partial {\mathit{\delta }}^{\mathrm{18}}\mathrm{O}}{\partial T_{\mathrm{c}}}$, does not change dramatically, it is clearly variable, as is $\frac{\partial {\mathit{\delta }}^{\mathrm{18}}\mathrm{O}}{\partial T_{\mathrm{0}}}$. The slopes $\frac{\partial d_{\mathrm{xs}}}{\partial T_{\mathrm{c}}}$ and $\frac{\partial d_{\mathrm{xs}}}{\partial T_{\mathrm{0}}}$ (comparable to β₁ and β₂ in Eq. 3, respectively) are highly variable. Indeed, the d_xs surface in Fig. 3 has a saddle at moderate condensation temperatures, over which $\frac{\partial d_{\mathrm{xs}}}{\partial T_{\mathrm{c}}}$ changes sign (Fig. 4c). This shows that the assumption of constant β and γ parameters in isotope-based temperature reconstructions is incorrect and may be reasonable only under narrow ranges of ΔT_site and ΔT_source. For plausible changes in site temperatures, assuming a fixed γ, for example, may not only lead to errors in magnitude but even to errors in the sign of γ and ultimately ΔT. The use of a fixed γ can introduce spurious variability into temperature reconstructions, particularly of T₀, the evaporation source temperature.

Evidence for the critical change in the sign of $\frac{\partial d_{\mathrm{xs}}}{\partial T_{\mathrm{c}}}$ (Fig. 4c) can actually be seen directly in Antarctic 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 $\frac{\partial d_{\mathrm{xs}}}{\partial T_{\mathrm{c}}}$, not from differing source conditions between sites. This exposes a fundamental flaw in the assumptions used in traditional isotope temperature reconstructions.

Figure 6 Time series and cross-correlation matrices for eight different deep Antarctic ice-core sites. (a, b) δ¹⁸O; (c, d) d_xs; and (e, f) d_ln. In the time series plots, each record is colored by its Holocene average δ¹⁸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.

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

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

4 Nonlinear temperature reconstructions

4.1 Reconstruction method

For every pair of T₀ and T_c inputs to SWIM there is a corresponding modeled value of δ¹⁸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., δ¹⁸O and d_ln. This defines a set of maps, with x and y axes of δ¹⁸O and d_ln and with z axes T_c and T₀, as shown in Fig. 7. To reconstruct T_c and T₀, the inverted model results may be used as a lookup table: a pair of δ¹⁸O and d_ln measurements determine a pair of T_c and T₀ 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.

Figure 7 Inverted T₀ and T_c surfaces as a function of modeled d_ln and δ¹⁸O of precipitation. (a) Surface shading and contours show condensation site temperature, T_c, in degrees Celsius (^∘C) as a function of d_ln and δ¹⁸O of precipitation at a site. (b) Surface shading and contours show evaporation source temperature, T₀, in degrees Celsius (^∘C) as a function of d_ln and δ¹⁸O of precipitation at a site.

The boundary conditions T_c and T₀ 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 δ¹⁸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, δ¹⁸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₀) instead of source region air temperature (T₀), since we assume climatological relationships between them. Indeed, RH₀ is a strong lever on the deuterium excess of evaporation. We reconstruct T₀ 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₀, and RH₀, by modeling multidimensional parameter spaces. This of course requires additional measured constraints, which could include the δ¹⁷O of precipitation, the accumulation rate, or other variables modeled in this framework, and is discussed in the Appendix.

4.2 Absolute temperature reconstructions

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

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

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

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

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

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

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

4.3 Comparison of linear and nonlinear reconstruction techniques

Linear temperature reconstruction using SWIM

We evaluate the significance of our approach by comparing our nonlinear reconstructions of condensation and evaporation temperature with reconstructions following the traditional linear approach. We calculate the equivalent linear β and γ coefficients for each of eight ice-core sites by regressing the T_c and T₀ temperature fields to subsets of δ¹⁸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, β₂ changes substantially between the glacial and the Holocene.

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

In the case of the surface temperature reconstruction, the errors introduced by linearization can be up to ±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₀ (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).

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

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

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

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

5 Discussion

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

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

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

We create a stack of each reconstructed temperature variable (the evaporation temperature T₀, 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.

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₀ from the ice-core reconstructions are shown as colored circles as in (a) for reference.

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

In Fig. 14 we plot the magnitude of warming since the LGM of Antarctic moisture source evaporation air temperatures for all ice-core records as a function of the mean latitude of the moisture source distribution for each site (based on a water-tagged general circulation model (GCM) simulations; see the Appendix). Additionally, we calculate the change in sea surface temperature (SST) during evaporation (red circles Fig. 14a) using the T₀–SST relationship from our model for comparison to other SST proxy reconstructions. While plotted as points, note that these changes in moisture source temperature reflect the integrated warming over the moisture source area. Modern moisture source distributions for each site are indicated by relative histograms along the latitude axis in Fig. 14a (see Fig. 1 for further information), which are colored to correspond to the ice-core site. The moisture source points are plotted at the longitude of the respective ice-core sites in Fig. 14b, though in reality the moisture source distributions (MSDs) have a significant meridional extent (see extended data in Fig. 8, Buizert et al., 2018) that is often asymmetrically to the west owing to the westerly winds. Changes in T₀ 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₀ reconstructions have low spatial resolution owing to broad moisture source distributions, they benefit from the temporal resolution and precision of the ice-core age scales compared to other proxy records of temperature from the Southern Hemisphere midlatitudes.

It is clear from Fig. 14 that the Antarctic continent warmed 2 to 3 times as much as the Southern Hemisphere midlatitude moisture source areas since the Last Glacial Maximum. This result is in line with other paleoclimate reconstructions, as well as modeling of the pattern of polar amplification since the LGM (Masson-Delmotte et al., 2006; Otto-Bliesner et al., 2006). In particular, our estimates of moisture source region changes agree with completely independent estimates from the MARGO compilation of SST changes (MARGO Project 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 reconstructions show more warming in the ice-core moisture source areas equatorward of the polar front than the PMIP3 mean and less warming over the surface of West Antarctica. We note that the magnitude and pattern of modeled Antarctic surface warming is predominately a function of imposed changes in ice sheet surface elevation to PMIP3 experiments. The extreme warming seen in parts of West Antarctica in the PMIP3 model mean (e.g., >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 ice-core records. Understanding the full set of processes responsible for the reconstructed pattern of Southern Hemisphere polar amplification is a topic for future work.

6 Conclusions

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

Qualitatively, δ¹⁸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 δ¹⁸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-isotope–temperature relationship (e.g., Jouzel et al., 1997). These discussions are conceptually useful. However, while space and time are obvious coordinates through which to understand climate, they are not the most relevant for water-isotope fractionation. Neither space nor time can independently cause water to change phase and fractionate.

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

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

Appendix A: Simple Water Isotope Model

The Simple Water Isotope Model (SWIM) is based on existing numerical Rayleigh-type distillation models (Merlivat and 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₀, 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₀ and RH₀ given a specified initial air temperature, T₀, 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₀, our model uses values of SST₀ and RH₀ based on fits to the modern climatology. We use three methods to calculate the climatological relationships over the interval −10^∘ $\le T_{\mathrm{a}}\le$ 28 ^∘C: a cubic spline with specified noise tolerance, the mean and median of SST and RH distributions within binned values of T_a, and 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₀ for different seasons and somewhat larger seasonal changes in the T_a and RH₀ relationship. We use the annual average fits hereafter and test the sensitivity of the water-isotope values of evaporation to these seasonal differences.

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.

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:

$$\begin{array}{}\text{(A1)}& {\mathrm{RH}}_n={\displaystyle \frac{\mathrm{RH}\times e_{\mathrm{s}}\left(T_{\mathrm{a}}\right)}{e_{\mathrm{s}}\left(\mathrm{SST}\right)}},\end{array}$$

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₀, and specified surface pressure, P₀, 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,

$$\begin{array}{}\text{(A2)}& r_{\mathrm{s}}={\displaystyle \frac{R_{\mathrm{d}}}{R_{\mathrm{wv}}}}\times {\displaystyle \frac{e_{\mathrm{s}}}{P-e_{\mathrm{s}}}},\end{array}$$

as functions of T, where $\frac{R_{\mathrm{d}}}{R_{\mathrm{wv}}}$ is the ratio of gas constants of dry air and water vapor.

We consider air parcels with mixed ice and liquid condensate (Ciais and Jouzel, 1994), in which the ice fraction smoothly increases as temperatures decreases below freezing. Many models, including 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 temperature-dependent functions for the cloud ice fraction derived from satellite observations (Hu et al., 2010) over the Southern Ocean and over land snow and ice surfaces, which preserve significantly more liquid water at colder temperatures than previous parameterizations (e.g., Kavanaugh and 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_{\mathrm{s}\left(\mathrm{eff}\right)}=r_{\mathrm{s}\left(\mathrm{ice}\right)}{F}_{\left(\mathrm{ice}\right)}+r_{\mathrm{s}\left(\mathrm{liq}\right)}{F}_{\left(\mathrm{liq}\right)}$, where r_s(ice) and r_s(liq) are the saturated mixing ratios of ice and liquid, respectively, and F_(ice) and F_(liq) are the (temperature-dependent) fractions of each phase of condensate.

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

At very cold temperatures moisture is removed not at saturation but at a specified level of supersaturation. The presence of both ice and liquid condensate in the cloud dictates a supersaturation of vapor over ice due to the difference in liquid and ice vapor pressures (Jouzel and 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 $\frac{{}^{\mathrm{2}}\mathrm{H}}{{}^{\mathrm{1}}\mathrm{H}}$, $\frac{{}^{\mathrm{18}}\mathrm{O}}{{}^{\mathrm{16}}\mathrm{O}}$, and $\frac{{}^{\mathrm{17}}\mathrm{O}}{{}^{\mathrm{16}}\mathrm{O}}$ ratios in water. We use conventional notation in which R is the number ratio of heavy to light isotopes of a species, for example ${}^{D}R=\frac{{}^{\mathrm{2}}\mathrm{H}}{{}^{\mathrm{1}}\mathrm{H}}$ and ${}^{\mathrm{18}}R=\frac{{}^{\mathrm{18}}\mathrm{O}}{{}^{\mathrm{16}}\mathrm{O}}$.

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

$$\begin{array}{}\text{(A3)}& {}^{\mathrm{18}}{\mathit{\alpha }}_{l-v}={\displaystyle \frac{{\left(\frac{{}^{\mathrm{18}}\mathrm{O}}{{}^{\mathrm{16}}\mathrm{O}}\right)}_{\mathrm{liquid}}}{{\left(\frac{{}^{\mathrm{18}}\mathrm{O}}{{}^{\mathrm{16}}\mathrm{O}}\right)}_{\mathrm{vapor}}}}={\displaystyle \frac{{}^{\mathrm{18}}{R}_{\mathrm{l}}}{{}^{\mathrm{18}}{R}_{\mathrm{v}}}}.\end{array}$$

We use the empirically determined, temperature-dependent equilibrium fractionation factors between liquid and vapor, ¹⁸α_eq(l−v) and ^Dα_eq(l−v), as well as those between vapor and ice, ¹⁸α_eq(i−v) and ^Dα_eq(i−v) (Majoube, 1970, 1971; Merlivat and Nief, 1967; Criss, 1999), with updates for the ice–vapor equilibrium fractionation factor found by Lamb et al. (2017).

A2.1 Evaporation from the ocean

The isotopic values of vapor evaporating from the ocean are determined, in part, by the isotopic values of the seawater. By definition, globally averaged seawater has δ values near 0 ‰. However, the δ¹⁸O of seawater (δ¹⁸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 δ¹⁸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 δ¹⁸O_sw (which changes with mean climate) and investigate the sensitivity of the model to these initial conditions. We assume a ¹⁷O_xs of seawater equal to 0, where ${}^{\mathrm{17}}{\mathrm{O}}_{\mathrm{xs}}=\mathit{\delta }{{}^{\prime }}^{\mathrm{17}}\mathrm{O}-\mathrm{0.528}\times \mathit{\delta }{{}^{\prime }}^{\mathrm{18}}$O, which is typically reported in per meg.

The atmosphere above the global oceans is not at saturation on average, with relative humidity typically around 80 % (Hartmann, 2015). Because of this 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

$$\begin{array}{}\text{(A4)}& {\mathit{\alpha }}_{\mathrm{diff}}={\left({\displaystyle \frac{D}{D^{*}}}\right)}^n,\end{array}$$

where D and D^* are the diffusivities of the light and heavy isotopes, respectively (Merlivat and Jouzel, 1979). The exponent n ranges from 0 to 1 and relates to the wind regime and speed and the ratio of turbulent to molecular diffusion. For the diffusive fractionation between ${\mathrm{H}}_{\mathrm{2}}^{\mathrm{18}}\mathrm{O}$ and ${\mathrm{H}}_{\mathrm{2}}^{\mathrm{16}}\mathrm{O}$ during initial evaporation, the fractionation factor ¹⁸α_diff equals 1.0 for pure turbulence and 1.0028 for pure molecular diffusion (Merlivat and Jouzel, 1979; Barkan and Luz, 2007).

Following Kavanaugh and Cuffey (2003), we do not explicitly consider surface wind speeds. Instead we use an empirical approach and the results of Uemura et al. (2008, 2010) for ¹⁸α_diff, who estimate the parameter based on measurements of δD, δ¹⁸O, and δ¹⁷O in vapor above the Southern Ocean. Uemura et al. (2010) find a value of ${}^{\mathrm{18}}{\mathit{\alpha }}_{\mathrm{diff}}=\mathrm{1.007}\pm \mathrm{0.0013}$ and 1.008±0.0018 when optimizing for observations of d_xs and ¹⁷O_xs of vapor, respectively. These results are within uncertainty of each other and of independent analysis by Pfahl and Wernli (2009), which found a value of 1.0076. Using a compilation of vapor measurements (including Uemura et al., 2008, 2010; Liu et al., 2014; Kurita et al., 2016; Benetti et al., 2017), we find that ¹⁸α_diff=1.009 leads to a good match between modeled and observed values of both d_xs and ¹⁷O_xs of vapor when using the observed values of T₀, 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).

$$\begin{array}{}\text{(A5)}& {\mathit{\varphi }}_{\mathrm{diff}}={\displaystyle \frac{{}^D{\mathit{\alpha }}_{\mathrm{diff}}-\mathrm{1}}{{}^{\mathrm{18}}{\mathit{\alpha }}_{\mathrm{diff}}-\mathrm{1}}}\end{array}$$

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 ¹⁸α_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 $\frac{{}^{\mathrm{17}}\mathrm{O}}{{}^{\mathrm{16}}\mathrm{O}}$ we use the following relationships, which are backed by both theory and empirical observation (Barkan and Luz, 2005, 2007): ${}^{\mathrm{17}}{\mathit{\alpha }}_{\mathrm{eq}}{=}^{\mathrm{18}}{\mathit{\alpha }}_{\mathrm{eq}}^{\mathrm{0.529}}$ for vapor and liquid in equilibrium and ${}^{\mathrm{17}}{\mathit{\alpha }}_{\mathrm{diff}}{=}^{\mathrm{18}}{\mathit{\alpha }}_{\mathrm{diff}}^{\mathrm{0.518}}$ for vapor diffusion.

An alternative to the empirical approach for diffusive fractionation outlined above is to use new results from kinetic molecular theory (Hellmann and Harvey, 2020), which lead to physically constrained temperature dependence of the diffusivity ratios for the isotopologues. These results lead to a temperature dependence of ϕ_diff that is different than the experimental results of Luz et al. (2009), being less variable and closer to the 0.88 value of Merlivat (1978). Using this formulation in the evaporation scheme requires us to choose a representative value of n (see Eq. A4) for the wind turbulence regime to match observations. A value of n between 0.22 and 0.32 leads to a range of results equivalent to that resulting from ¹⁸α_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),

$$\begin{array}{}\text{(A6)}& {\mathit{\alpha }}_{\mathrm{evap}}={\displaystyle \frac{R_{\mathrm{o}}}{R_{\mathrm{e}}}}={\displaystyle \frac{{\mathit{\alpha }}_{\mathrm{eq}}{\mathit{\alpha }}_{\mathrm{diff}}(\mathrm{1}-{\mathrm{RH}}_n)}{\mathrm{1}-{\mathit{\alpha }}_{\mathrm{eq}}{\mathrm{RH}}_n\left(\frac{R_{\mathrm{v}}}{R_{\mathrm{o}}}\right)}},\end{array}$$

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):

$$\begin{array}{}\text{(A7)}& R_{\mathrm{v}}={\displaystyle \frac{R_{\mathrm{o}}}{{\mathit{\alpha }}_{\mathrm{eq}}\times \left({\mathit{\alpha }}_{\mathrm{diff}}+{\mathrm{RH}}_n\left(\mathrm{1}-{\mathit{\alpha }}_{\mathrm{diff}}\right)\right)}}.\end{array}$$

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

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

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

We investigate the influence of the closure assumption in a few ways. First, we examine closure globally rather than locally. The mean ocean has δ values of about 0 ‰, and global average precipitation has δ¹⁸O = −4.5 ‰ and δD $=-\mathrm{26.7}$ ‰ (Craig and Gordon, 1965). Considering the global average steady state in which the ocean is the ultimate vapor source for precipitation (Merlivat and Jouzel, 1979; Criss, 1999) , the delta values of precipitation must reflect the net loss by evaporate from the ocean. Thus, globally, the ratio $\frac{R_{\mathrm{o}}}{R_{\mathrm{e}}}$ is 1.0267 for D and 1.0045 for ¹⁸O. Instead of equating R_e and R_v locally, we define a global ${\stackrel{\mathrm{‾}}{\mathit{\alpha }}}_{\mathrm{evap}}={\left(\frac{R_{\mathrm{o}}}{R_{\mathrm{e}}}\right)}_{\mathrm{global}}$ and solve for R_v. Substitution into Eq. (A6) and rearranging leads to

$$\begin{array}{}\text{(A8)}& R_{\mathrm{v}}=R_{\mathrm{o}}\left(\mathrm{1}-{\displaystyle \frac{{\mathit{\alpha }}_{\mathrm{eq}}{\mathit{\alpha }}_{\mathrm{diff}}\left(\mathrm{1}-{\mathrm{RH}}_n\right)}{{\stackrel{\mathrm{‾}}{\mathit{\alpha }}}_{\mathrm{evap}}}}\right){\left({\mathit{\alpha }}_{\mathrm{eq}}{\mathrm{RH}}_n\right)}^{-\mathrm{1}}.\end{array}$$

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

Next, we consider specific mixing of evaporative conditions instead of the generalized globally mixed case above. Rather than the local ocean being the only source of vapor, we can consider a simple scenario in which the isotopic composition of the boundary layer, R_v, is comprised of both local evaporate, R_e, and vapor evaporated at some distal location and advected to the site, $R_{\mathrm{v}}=(\mathrm{1}-\mathit{\theta })R_{\mathrm{e}}+\mathit{\theta }{R}_{\mathrm{distal}}$, where θ is the fraction of nonlocal vapor with composition R_distal. The local evaporative conditions are defined by T₀, while the distal conditions may be either warmer or colder than T₀. 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₀, with a full range of contributions from both a 5 ^∘C warmer moisture source (T₀+5 ^∘C, red) and a 5 ^∘C colder moisture source (T₀−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.

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

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

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

It is important to note that in Fig. A2 we use climatological correlations between T₀, SST, and RH, while the observational data represent far shorter time intervals, mostly from one season. When using the observed values of T₀, SST, and RH at the time of the observational measurements (Uemura et al., 2008; Liu et al., 2014), we are able to capture the complex variability in the isotopic values of the vapor on those given days. For example, in Fig. A4 we compare modeled d_xs of vapor and modeled ¹⁷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 ¹⁷O_xs with SST₀ and RH₀ 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₀, SST₀, and RH₀, 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 δ¹⁸O_sw of the ocean (Fig. A5g–i) as well as the value of α_diff (Fig. A5j–l).

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

Figure A6 Relationship between observed and model d_xs (a) and ¹⁷O_xs (b) of vapor over the Southern Ocean. Observed vapor values are from Uemura et al. (2008, 2010) (large dots) and Liu et al. (2014) (small dots). Modeled values use the reported T₀, SST₀, and RH₀ from the observations and four different values of ¹⁸α_diff (shown in color of dots). The 1:1 line is shown in black.

The direct comparisons of observed and modeled vapor composition using observed T₀, SST₀, and RH₀ 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 temporal variability in δ¹⁸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

$$\begin{array}{}\text{(A9)}& {\displaystyle \frac{d\ln \left(R\right)}{d\ln \left(f\right)}}=\mathit{\alpha }-\mathrm{1}\end{array}$$

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,

$$\begin{array}{}\text{(A10)}& f={\displaystyle \frac{q}{q_{\mathrm{0}}}}={\displaystyle \frac{r_{\mathrm{s}}}{r_{\mathrm{s}\mathrm{0}}}}.\end{array}$$

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

$$\begin{array}{}\text{(A11)}& d\ln \left(R\right)=({\mathit{\alpha }}_{\mathrm{tot}}-\mathrm{1})d\ln \left(f\right).\end{array}$$

The kinetic fractionation factor, α_k, depends on the supersaturation of vapor over ice, S_i:

$$\begin{array}{}\text{(A12)}& {\mathit{\alpha }}_{\mathrm{k}}={\displaystyle \frac{S_{\mathrm{i}}}{{\mathit{\alpha }}_{\mathrm{eq}}\times \frac{D}{D^{*}}\left(S_{\mathrm{i}}-\mathrm{1}\right)+\mathrm{1}}}.\end{array}$$

Following Jouzel and Merlivat (1984), we use the ratio of diffusivities for oxygen isotopes $\frac{D^{\mathrm{16}}}{D^{\mathrm{18}}}=\mathrm{1.0285}$ during moisture transport, representative of pure molecular diffusion and ignoring the negligible ventilation effect. Likewise we use $\frac{D^{\mathrm{1}}}{D^{\mathrm{2}}}=\mathrm{1.0251}$ for the ratio of diffusivities of hydrogen isotopes. These values imply a constant ϕ_diff during transport equal to 0.88 (Jouzel and Merlivat, 1984) rather than the 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),

$$\begin{array}{}\text{(A13)}& {\mathit{\alpha }}_{\mathrm{eff}}={\mathit{\alpha }}_{\mathrm{tot}(\mathrm{l}-\mathrm{v})}{F}_{\left(\mathrm{liq}\right)}+{\mathit{\alpha }}_{\mathrm{tot}(\mathrm{i}-\mathrm{v})}{F}_{\left(\mathrm{ice}\right)}.\end{array}$$

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:

$$\begin{array}{}\text{(A14)}& S_{\mathrm{i}}=a-b\times T,\end{array}$$

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 $\frac{\mathrm{d}\ln \left(f\right)}{\mathrm{d}T}$) and the supersaturation driving kinetic fractionation relates temperature to δ values (largely through Eq. A12), an inconsistency in the model's view of supersaturation can influence the modeled water-isotope–temperature relationship in unphysical ways.

To resolve this physical inconsistency, precipitation only occurs in SWIM when the parcel reaches the prescribed supersaturation by dictating an effective saturated mixing ratio of the air parcel, in which 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 δ¹⁸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 $\frac{\partial {\mathit{\delta }}^{\mathrm{18}}\mathrm{O}}{\partial \mathit{\delta }\mathrm{D}}$ that is not compatible with measured values in Antarctic precipitation. This gives additional credence to sustained supersaturation at cold temperatures.

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

A3 Application to Antarctic ice-core sites

A3.1 Moisture source distributions

The atmospheric circulation transports moisture poleward of ≈30^∘ S (Fig. 1a). The mean evaporation latitude of moisture that precipitates at any given site can be estimated from moisture-tagged GCM experiments (Markle et al., 2017). The difference between the mean latitudes of moisture evaporation and precipitation steadily increases between the subtropics and the pole (Fig. 1a and b). The mean evaporative latitude of moisture that precipitates in Antarctica ranges between 44 and 50^∘ S (Fig. 1c). The surface elevation of the ice sheet is a strong predictor of the mean latitude of precipitation, with higher-elevation sites having more equatorward moisture sources (Fig. 1d) due to topographic isolation (Sodemann and 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 relationship between the condensation temperature and the surface temperature (T_s) over Antarctica and compare the surface temperature to the weighted annual mean condensation temperature in both ERA-40 reanalysis (1980–2002) and MAR, a high-resolution mesoscale model forced by ERA-40 (see Fig. 8, Masson-Delmotte et al., 2008). In both models the upper bound of the Antarctic condensation temperature appears to be set by the peak inversion temperature, though condensation temperatures are on average colder than the peak inversion temperature (meaning simply that condensation occurs at a range of temperature up to the peak inversion temperature). Masson-Delmotte et al. (2008) calculate a best fit of the surface to condensation temperature slope of 0.65 ^∘C ^∘C⁻¹ in the ERA-40 data, consistent with previous work that found a slope of 0.67 ^∘C ^∘C⁻¹ (Connolley, 1996; Jouzel and Merlivat, 1984). The spread of condensation temperatures in the higher-resolution MAR model suggests colder condensation temperatures than in the lower-resolution reanalysis (Masson-Delmotte et al., 2008). The strength of the Antarctic inversion diminishes with increasing surface temperature (Connolley, 1996), and relatively warm Antarctic surface temperatures (e.g., $>-\mathrm{20}$ ^∘C) are associated with condensation temperatures colder than the surface temperature (Masson-Delmotte et al., 2008).

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

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

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

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

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

Examining all modeled condensation temperatures for samples in the MD08 and GNIP datasets with reported surface temperatures below 15 ^∘C, we find slopes between 0.62 and 0.67 ^∘C ^∘C⁻¹. 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⁻¹ (Fig. A9), depending on the model assumptions (0.69 ^∘C ^∘C⁻¹ 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⁻¹, depending on whether we average samples from the last 50 or 100 years, though only five points describe these lines.

Based on the above results we use the equation $T_{\mathrm{c}}=\mathrm{0.69}\frac{{}^{\circ }\mathrm{C}}{{}^{\circ }\mathrm{C}}{T}_{\mathrm{s}}-\mathrm{8.2}{}^{\circ }\mathrm{C}$ as our base estimate to reconstruct Antarctic surface temperatures,; however, we consider an uncertainty of $\pm \mathrm{0.02}\frac{{}^{\circ }\mathrm{C}}{{}^{\circ }\mathrm{C}}$ in the slope. Our base estimate leads to good agreement with the observed relationship between global δ¹⁸O and surface temperature (Fig. A10). We show our temperature-dependent isotope slopes in Fig. A11.

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

A3.3 Seasonality

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

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

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

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

Figure A11 Partial derivatives of isotope state spaces with respect to surface temperature: (a) $\frac{\partial {\mathit{\delta }}^{\mathrm{18}}\mathrm{O}}{\partial T_{\mathrm{s}}}$, (b) $\frac{\partial \mathit{\delta }\mathrm{D}}{\partial T_{\mathrm{s}}}$, (c) $\frac{\partial d_{\mathrm{xs}}}{\partial T_{\mathrm{s}}}$, (d) $\frac{\partial d_{\ln }}{\partial T_{\mathrm{s}}}$. Shading and contours in all panels is the slope in ‰ ^∘C⁻¹.

A4 Tuning the Simple Water Isotope Model

We tune the Simple Water Isotope Model by adjusting the temperature dependence of the supersaturation of vapor over ice. Given insufficient observational and physical constraints, we parameterize the supersaturation as a linear function of temperature (Jouzel and Merlivat, 1984) as above, $S_{\mathrm{i}}=a+b\times T$, set a=1, and tune the slope, b. The supersaturation has a strong influence on the kinetic fractionation (Eq. A12) and thus the relationship between δD and δ¹⁸O in vapor and precipitation. We tune the model to yield the observed relationship between δD and δ¹⁸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 ¹⁷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). Data are reported relative to the VSMOW (Vienna Standard Mean Ocean Water) standard and normalized to SLAP.

The global relationship between δD and δ¹⁸O has an approximate slope of 8, as codified in the historical definition of the deuterium excess parameter (Dansgaard, 1964). However, the slope is not fundamental (Craig, 1961); as discussed in Sect. 1.2 the true observed relationship is nonlinear (Uemura et al., 2012), as is the theoretical relationship even in the absence of kinetic fractionation (Markle et al., 2017). Uemura et al. (2012) find an empirical fit between δ^′D and $\mathit{\delta }{{}^{\prime }}^{\mathrm{18}}\mathrm{O}$ in a global dataset for precipitation. They use a 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 δ¹⁸O, given any amount of distillation, depends on the ratio of $\frac{\text{exp}{(}^D{\mathit{\alpha }}_{\mathrm{tot}})}{\text{exp}{(}^{\mathrm{18}}{\mathit{\alpha }}_{\mathrm{tot}})}$, where each α_tot is itself a nonlinear function of temperature as outlined throughout Sect. A2.2. The ratio of exponential functions can be estimated to any arbitrary degree of accuracy with a polynomial function.

Our modern dataset includes several new sets of measurements in addition to those used in Uemura et al. (2012). We find similar coefficients in a second-order polynomial fit between $\mathit{\delta }{{}^{\prime }}^{\mathrm{18}}$O and δ^′D in our larger dataset compared to those found by Uemura et al. (2012): $A=-\mathrm{29.2}$ and B=8.45 (see Eq. 4). Because these coefficients are not significantly different than those previously published and for consistency with that work, we use the coefficients found by Uemura et al. (2012) ($A=-\mathrm{28.5}$ and B=8.47) to define a logarithmic deuterium excess parameter, d_ln. We find no benefit or justification for using 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₀ varies from 0 to 28 ^∘C; T_c from 27 to −60 ^∘C). We then compare the resulting cloud of modeled $\mathit{\delta }{{}^{\prime }}^{\mathrm{18}}\mathrm{O}$, δ^′D, and d_ln, finding a second-order polynomial fit between the modeled δ^′D and $\mathit{\delta }{{}^{\prime }}^{\mathrm{18}}\mathrm{O}$ from the ensemble of temperature trajectories. We tune the model by iteratively adjusting the b value in the supersaturation parameterization to minimize the difference between the modeled and observed relationship between δ^′D and $\mathit{\delta }{{}^{\prime }}^{\mathrm{18}}\mathrm{O}$ (Eq. 4). This is easily visualized in a plot of $\mathit{\delta }{{}^{\prime }}^{\mathrm{18}}\mathrm{O}$ and d_ln (e.g., Fig. A12a), as the average $\mathit{\delta }{{}^{\prime }}^{\mathrm{18}}\mathrm{O}$-to-d_ln relationship is flat in measured samples by definition.

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

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

The appropriate weighting of model realizations with different T₀ 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₀, though this is not a strong constraint on the tuning. While one can reasonably reject most possible values of b (as in Fig. A13), we cannot justifiably discern between others within a small range (e.g., $\mathrm{0.005}\le b\le \mathrm{0.0055}$). This gives a small but inherent uncertainty to the model tuning and in turn our temperature reconstructions. In spite of these limitations, the tuning procedure reproduces the observed isotope relationships well.

Figure A12 Tuning SWIM with linear supersaturation parameterization. (a) The modeled δ¹⁸O and d_ln of precipitation (colored circles) for a range of condensation and evaporation temperatures. Color shading shows source region evaporation temperature in degrees Celsius (^∘C). Black dots are the target dataset as described in the text. Modeled results are for the optimized supersaturation parameterization ($S_{\mathrm{i}}=\mathrm{1}-b\times T$, b=0.00525 ^∘C⁻¹) using the local closure assumption and the NCEP/NCAR reanalysis dataset for source region correlations. (b) Same as panel (a) for the modeled and target δ¹⁸O and δD of precipitation. (c) Same as panel (a) for the modeled and target δ¹⁸O and d_xs of precipitation.

Figure A13 Comparison of tuned and rejected supersaturation parameterizations, $S_{\mathrm{i}}=\mathrm{1}-b\times T$. (a) Tuned example: b=0.00525 ^∘C⁻¹. Same as in Fig. A12a. The modeled δ¹⁸O and d_ln of precipitation (colored circles) and target dataset (black dots). Color shading shows source region evaporation temperature in degrees Celsius (^∘C). (b) Same as panel (a), but for a rejected tuning: b=0.003 ^∘C⁻¹. The modeled d_ln values curve upward strongly with δ¹⁸O, incongruently with the target data. (c) Same as panel (a), but for another rejected tuning: b=0.007 ^∘C⁻¹. The modeled d_ln values curve downward strongly with δ¹⁸O and incongruently with the target data.

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

While the linear, temperature-dependent parameterization of supersaturation is both simple and widely used, the physical processes determining supersaturation are complex. To understand the sensitivity of the model to this parameterization we also test a nonlinear parameterization of supersaturation, $S_{\mathrm{i}}=a-b\times T-c\times T^{\mathrm{2}}$. If the physical source of the high supersaturation at very cold conditions is related to the absence of condensation nuclei, the supersaturation may not linearly increase at very cold temperatures, as there could be diminishing returns as the atmosphere becomes increasingly clean. A small but positive c parameter that gradually decreases the slope in S_i with decreasing temperature could be plausible. Only very small values of the second-order term, c, those of order $\mathrm{5}\times {\mathrm{10}}^{-\mathrm{6}}$ ^∘C⁻², are reconcilable with the observed $\mathit{\delta }{{}^{\prime }}^{\mathrm{18}}\mathrm{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 moisture-weighted average of the unmixed parcels.

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

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

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

We use the simple model to assess the range of final isotopic values of precipitation that can arise from two air parcels, which evaporated at the same initial conditions, then mixed at a range of different temperatures during transport. Air parcels are evaporated at specified initial air temperature (T₀=10 ^∘C), cooled, and randomly mixed at any combination of temperatures whose difference does not exceed a threshold, then cooled the remainder of the temperature pathway to −30 ^∘C. We assume no preferential temperature of mixing and no preferential difference in temperature during mixing, though we assume a normally distributed probability of mixing with temperature differences up to 5 ^∘C (a high-end estimate based on the above analysis of daily temperature differences). This process is repeated 10 000 times to create a distribution of final δ¹⁸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₀ 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 δ¹⁸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 δ¹⁸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.

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

A6 Optimal coordinates for reconstruction technique

Consider a water sample with mean values of δ¹⁸O and δD and normally distributed uncertainties σ₁₈ 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 σ₁₈ and σ_D. How do these uncertainties influence the temperature reconstruction?

In Fig. A15 we examine the estimation of T₀ using coordinates of δ¹⁸O and δD (Fig. A15a), δ¹⁸O and d_xs (Fig. A15b), and δ¹⁸O and d_ln (Fig. A15c). The uncertainties in the position of the measurement along both the x axis (δ¹⁸O) and y axis (either δD, d_xs, or d_ln) combine to give the total uncertainty in the position on the T₀ surface, shown as the targets in Fig. A15a–c. The total combined uncertainty in the estimation of T₀ is shown as probability density functions (PDFs) for each method in Fig. A15d. All estimates yield the same mean value of reconstructed T₀; 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₀ 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 δ¹⁸O values, that uncertainty will be spread across a wider range of isotherms. The axis of the influence of T₀ on d_xs (and δD) is rotated with respect to its axis of variability in d_xs.

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

This result is of course ultimately tied to the same reasons that 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 δ¹⁸O and d_ln measurements (corrected for changes in seawater δ¹⁸O) from the eight Antarctic deep ice-core records examined in this study.

Figure A16 Inverted T₀ and T_c surfaces as a function of modeled d_ln and δ¹⁸O of precipitation as in Fig. 7. Overlain are pairs of d_ln and δ¹⁸O measurements for eight different deep ice-core records, whose site locations are shown on the inset map. See Sect. 4.2 in the text for details on ice-core records.

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

Figure 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.

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.

A8 Correlation of nonlinear and linear reconstruction techniques

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

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

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

By accounting for the fundamental nonlinearities in water-isotope distillation we are able to reveal more coherent underlying climate signals in source region temperatures, which are otherwise obscured by linear temperature reconstruction techniques. For analogous reasons, Markle et al. (2017) argued that the logarithmic deuterium excess parameter 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₀ and d_ln parameters share the same correlation pattern amongst the ice cores and show substantially more coherence than either linearly reconstructed T₀ or d_xs. The correlation pattern of d_ln and d_xs between all core sites (Fig. A24) reveals how nonlinear effects alter the traditionally defined d_xs at the coldest Antarctic temperatures. The broad change from positive to negative correlation of d_ln to d_xs across sites is a reflection of the change in sign of $\frac{\partial d_{\mathrm{xs}}}{\partial T_{\mathrm{c}}}$ as a function of T_c.

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

Figure A22 The increase in R² (shared variance) of the nonlinear reconstructions of T₀ for all sites over the linear reconstructions (linearized around Holocene conditions). All records are ordered by their approximate modern surface temperature.

Figure A23 Correlation matrices for (a) δ¹⁸O, (b) δD, (c) d_xs, and (d) d_ln among all core sites. All records are ordered by their approximate modern surface temperature.

Figure A24 Correlation matrix between d_xs and d_ln between all core sites. All records are ordered by their approximate modern surface temperature.

A9 Temperature reconstruction uncertainty

In this section we investigate uncertainty in our temperature reconstructions by examining the sensitivity of our results to assumptions and parameterizations in the model. We can compare reconstructed T_c and T₀ from a set of δ¹⁸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₀ 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 δ¹⁸O and d_ln. Thus, uncertainty arising from a given parameter may vary between ice-core sites. In general, varying a parameter in the model results in patterns of the partial slopes in δ¹⁸O and d_ln with T_c and T₀ that are similar, but shifted in the T_c and T₀ space. A consequence of this is that the uncertainty in absolute values of reconstructed T_c and T₀ is generally larger than uncertainty in their relative variability.

Figure A25 Temperature reconstruction sensitivity to model parameterizations. T_c and T₀ reconstructions using WDC δ¹⁸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 ¹⁸α_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₀, constant RH along path), and all moisture is removed above fixed 80 % or 90 % RH.

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

A9.1 Sensitivity to model parameters

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

Aspects of the initial evaporation scheme introduce uncertainty into our reconstructions. The value of ${\mathit{\alpha }}_{\mathrm{diff}}^{\mathrm{18}}$ during evaporation is important in setting the initial isotopic values of vapor. While we find the value 1.009 to give the best fit to modern observations, values within a small range may be defensible (Fig. A6). The local closure assumption used in the evaporation scheme has known limitations (Risi et al., 2010), representing an 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₀ up to 1.5 ^∘C for the WDC record, while differences in absolute T_c are smaller (≤1 ^∘C). Relative variability in T₀ and T_c is similar when using either closure assumption and ≤0.3 ^∘C.

We also examine the influence of the source relative humidity parameterization on our temperature reconstructions. In our base model we use climatological correlations to determine an initial relative humidity given an initial air temperature; colder surface air temperatures over the ocean are associated with slightly higher relative humidity. We show the difference in reconstructions due to using either the NCEP or ERA-Interim reanalysis. We can also ask how different our reconstructed T_c and T₀ in WDC would be if we used fixed mean values of initial relative humidity rather than values that depend on T₀. 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₀ and T_c as functions of δ¹⁸O and d_ln to maps produced by a large ensemble of model runs that incorporate stochastic mixing. We consider a range of final condensation temperatures from −20 to −50 ^∘C. We generate random pairs of pseudo-adiabatic cooling pathways ending at every value of T_c and random values of T₀ 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⁴ unmixed and 7.5×10⁵ 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₀ and T_c as functions of δ¹⁸O and d_ln (at a resolution of 0.1 ‰, Fig. A26). Due to the stochastic mixing these maps are unevenly populated. The potential influence of mixing on our reconstruction technique can be seen in the difference between the mean T_c and T₀ maps resulting from the mixed (denoted with the subscript M) and unmixed pathways (Fig. A26). In Fig. A26 we show the histograms of $\mathrm{\Delta }{T}_{\mathrm{c}}=T_{\mathrm{c},\mathrm{M}}-T_{\mathrm{c}}$ and $\mathrm{\Delta }{T}_{\mathrm{0}}=T_{\mathrm{0},\mathrm{M}}-T_{\mathrm{0}}$. We test a range of mixing and threshold values in multiple Monte Carlo simulations. In all cases the mean values of $T_{\mathrm{c},\mathrm{M}}-T_{\mathrm{c}}$ and $T_{\mathrm{0},\mathrm{M}}-T_{\mathrm{0}}$ are very near zero ($<|\pm \mathrm{0.06}$ ^∘C| for T_c, and $<|\pm \mathrm{0.02}$ ^∘C| for T₀). 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₀ in all tests. We find similar results when using isobaric rather than pseudo-adiabatic cooling pathways. Including moisture sources with $T_{\mathrm{0}}<\mathrm{0}{}^{\circ }$C in our mixing analysis has no significant influence on the mean difference between the weighted mean maps of either T_c or T₀, though expanding the range of moisture sources to include T₀<0 ^∘C does increase the range of $T_{\mathrm{0},\mathrm{M}}-T_{\mathrm{0}}$ by over a degree, in agreement with the analysis of unmixed pathways.

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⁵ stochastically mixed distillation pathways, colored by T_c and T₀, respectively. (b, f) The water-isotope ratios of precipitation resulting from the corresponding 1.5×10⁴ unmixed distillation pathways, colored by T_c and T₀, respectively. (c, g) The difference in T_c and T₀, respectively, between the mixed and unmixed pathways as a function of the δ¹⁸O and d_ln space. (d, h) Histograms of the differences, ΔT_c and ΔT₀, from all points in the δ¹⁸O and d_ln space, resulting from mixing.

A9.3 Combined uncertainty estimates

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

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

Based on the tests in Sect. A2.1, a conservative estimate of the uncertainty arising from mixing at the evaporation source is half the absolute mean difference in reconstructions employing the local and global closure assumptions, applied symmetrically about the base scenario (local closure). Because of the stochastic nature of our atmospheric mixing simulations, our estimates of the T_c and T₀ differences are nonuniform and unevenly populated, making interpolation in the δ¹⁸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 δ¹⁸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_{\mathrm{c}}\propto \mathrm{0.69}\pm \mathrm{0.02}\frac{{}^{\circ }\mathrm{C}}{{}^{\circ }\mathrm{C}}{T}_{\mathrm{s}}$ to estimate this uncertainty, which is added in quadrature to the above uncertainties in T_s. An example of this spread of uncertainty in the WDC T_s reconstruction is shown in Fig. A27. Reconstructions of T_s, T_c, and T₀ for several major ice-core records along with the combined relative uncertainty are those reconstructions is shown in Fig. 9.

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_{\mathrm{c}}\propto \mathrm{0.69}\frac{{}^{\circ }\mathrm{C}}{{}^{\circ }\mathrm{C}}{T}_{\mathrm{s}}$), while the spread associated with a $\pm \mathrm{0.02}\frac{{}^{\circ }\mathrm{C}}{{}^{\circ }\mathrm{C}}$ uncertainty in the scaling factor is shown in red and gold. For comparison the results from a weaker slope of $\mathrm{0.65}\frac{{}^{\circ }\mathrm{C}}{{}^{\circ }\mathrm{C}}$ are shown in blue. This range of scaling factors has little impact on the reconstructed temperature history and may be difficult to see.

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₀ (^∘C) using base model assumptions. For initial evaporation air temperatures below 0 ^∘C there are non-unique results in the δ¹⁸O and d_ln space.

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 ^∘C), with the contribution from T₀<0 ^∘C weighted by the final r_s (red) and with a 90 % contribution from T₀>0 ^∘C as well as 10 % from T₀<0 ^∘C (gold). (b) Same as in (a) but for reconstructions of T₀.

A9.4 Uniqueness and source temperature

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

Those caveats aside, we investigate the sensitivity of our temperature reconstructions to this potential non-uniqueness in the water-isotope state spaces. We run SWIM through a large field of T₀ 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 ^∘C, 10 % T₀<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₀ variability: the standard deviation of differences in reconstructed T_c is less than 0.07 ^∘C using either method and less than 0.19 ^∘C for T₀. 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₀ 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 δ¹⁸O and d_ln space, as well as the actual values of δ¹⁸O and d_ln in ice-core measurements, the model must extrapolate to extremely cold T₀ values for the below 0 ^∘C side of the folded surface. These values of T₀ 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₀ when a much warmer one (and perhaps a reduced RH₀) is actually correct. The net result of these considerations is that the analysis above should represent a quite conservative estimate of the influence of non-uniqueness on our temperature reconstructions and their relative variability; the real effect is likely much smaller though difficult to quantify precisely.

A10 Comparison to previous reconstructions

We next reconstruct site and source temperatures for four East Antarctic ice-core records and compare to previously published linear reconstructions. We use records of δ¹⁸O and δD measurements (and calculate d_xs and d_ln) from the Vostok (Jouzel et al., 1997; Uemura et al., 2012), EPICA Dome C (EDC) (Stenni et al., 2004, 2010), EPICA Dronning Maud Land (EDML) (Stenni et al., 2010), and Dome Fuji records (Uemura et al., 2012). After seawater correction, we use the ice-core δD and d_xs for the linear reconstruction and δ¹⁸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₀, 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, $\frac{\partial d_{\mathrm{xs}}}{\partial T_{\mathrm{c}}}$ and $\frac{\partial d_{\mathrm{xs}}}{\partial T_{\mathrm{0}}}$. If one assumes these slopes are linear over a given range in T₀ 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 δ¹⁸O and the temperature boundary conditions, though the nonlinearities in these slopes are much smaller.

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

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

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

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

A11 Three-parameter reconstructions

In the approach outlined above, we have considered the boundary conditions T_c and T₀ to be the only independent input variables. In particular, we have assumed that the source region relative humidity, RH₀, is a dependent variable whose value is not fixed but determined by climatological correlation with T₀. 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₀ is dependent on T₀ and reconstruct three independent climate variables (T_c, T₀, and RH₀) if we have three independent constraints. While δ¹⁸O and d_ln alone are not sufficient, ${}^{\mathrm{17}}{\mathrm{O}}_{\mathrm{xs}}=\mathit{\delta }{{}^{\prime }}^{\mathrm{17}}\mathrm{O}-\mathrm{0.528}\mathit{\delta }{{}^{\prime }}^{\mathrm{18}}$O (Landais et al., 2008) can in principle provide the necessary additional information. We can allow T_c, T₀, and RH₀, to all vary as independent variables, defining a three-dimensional parameter space through which SWIM is run to produce three-dimensional isotope state spaces.

While promising, this method currently has practical limitations. Our model does not reproduce the observed ¹⁷O_xs-to-δ¹⁸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 ¹⁷O_xs in Antarctic precipitation may be too large to offer useful discrimination amongst variations in T_c, T₀, and RH₀ (Schoenemann et al., 2014).

The ¹⁷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 ¹⁷O_xs lead to large changes in the absolute value of reconstructed source region conditions (T₀ and RH₀). It is worth noting that absolute values of ¹⁷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 δ¹⁸O, d_ln, and ¹⁷O_xs, does not necessarily lead to unique values of the boundary conditions.

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