Snapshots of mean ocean temperature over the last 700 000 years using noble gases in the EPICA Dome C ice core

Together with the latent heat stored in glacial ice sheets, the ocean heat uptake carries the lion's share of glacial–interglacial changes in the planetary heat content, but little direct information on the global mean ocean temperature (MOT) is available to constrain the ocean temperature response to glacial–interglacial climate perturbations. Using ratios of noble gases and molecular nitrogen trapped in the Antarctic EPICA Dome C ice core, we are able to reconstruct MOT for peak glacial and interglacial conditions during the last 700 000 years and explore the differences between these extrema. To this end, we have to correct the noble gas ratios for gas transport effects in the firn column and gas loss fractionation processes of the samples after ice core retrieval using the full elemental matrix of N₂, Ar, Kr, and Xe in the ice and their individual isotopic ratios. The reconstructed MOT in peak glacials is consistently about 3.3 ± 0.4 ^∘C cooler compared to the Holocene. Lukewarm interglacials before the Mid-Brunhes Event 450 kyr ago are characterized by 1.6 ± 0.4 ^∘C lower MOT than the Holocene; thus, glacial–interglacial amplitudes were only about 50 % of those after the Mid-Brunhes Event, in line with the reduced radiative forcing by lower greenhouse gas concentrations and their Earth system feedbacks. Moreover, we find significantly increased MOTs at the onset of Marine Isotope Stage 5.5 and 9.3, which are coeval with CO₂ and CH₄ overshoots at that time. We link these CO₂ and CH₄ overshoots to a resumption of the Atlantic Meridional Overturning Circulation, which is also the starting point of the release of heat previously accumulated in the ocean during times of reduced overturning.

1 Introduction

Over the last million years, Earth's climate has experienced pronounced changes in global ice volume (Bintanja et al., 2005) and hence sea level, accompanied by significant temperature changes between cold glacials and warmer interglacials. These temperature changes are observed with different amplitudes both on land (EPICA community members, 2004; Tzedakis et al., 2006; Melles et al., 2012) and in the ocean (Elderfield et al., 2012; Shakun et al., 2012, 2015). In particular, due to the large size and heat capacity of the ocean reservoir, global mean ocean temperature (MOT) changes (together with the latent heat stored in waxing ice sheets) represent an integrated signal of Earth's energy imbalance in the past (Baggenstos et al., 2019). For example, today, the global ocean is taking up about 90 % of the excess heat from anthropogenic global warming, dominating the current changes in the planetary energy budget (Gleckler et al., 2016; von Schuckmann et al., 2020).

In contrast to today's warming, where radiative forcing is caused through anthropogenic emissions of CO₂, past climate cycles are generally believed to have been driven by orbital changes in the latitudinal and seasonal distribution of incoming solar radiation (Milankovic, 1941), with changing greenhouse gas concentrations representing an important amplifying feedback. Only recently did Tzedakis et al. (2017) propose a simple rule for the occurrence of full deglaciations based on the amount of effective incoming energy flux (defined as the caloric summer half-year insolation at 65^∘ N corrected by a term that takes into account the time elapsed since the previous deglaciation). While this rule is in agreement with the occurrence of all full deglaciations over the last million years, it does not intend to predict differences between individual interglacials in ice volume and climate. For example, Marine Isotope Stage (MIS) 11 is characterized by a very long interval of interglacial conditions and a sea level at least equally high as in MIS 5.5 and likely more than 6 m higher than in the Holocene (Dutton et al., 2015) . However, its effective energy is just reaching the threshold suggested by Tzedakis et al. (2017), while MIS 5.5 exhibits the largest effective energy over the last 800 kyr and the second largest over the last 2.6 Myr (Tzedakis et al., 2017). Moreover, lukewarm interglacials (i.e., interglacials significantly colder than the Holocene with increased ice volume) were observed in local temperature reconstructions and global sea level (EPICA community members, 2004; PAGES Past Interglacials Working Group, 2016; Shakun et al., 2015; Elderfield et al., 2012) in the time interval 800–450 kyr BP (before present; present is defined as 1950). Again, the effective energy during these lukewarm interglacials is not systematically lower than for the subsequent interglacials, requiring additional energetic or radiative changes in the climate system to explain the difference between the individual interglacials. MOT may play an important role in explaining these features, as a higher ocean heat content could represent a compensating mechanism for the excess energy not used for the melting of the larger remnant ice sheets during lukewarm interglacials. The opposite is also true: if lower MOT parallels larger ice volume during lukewarm interglacials, changes in the radiative balance of the Earth are required to explain the differences between the interglacials before and after 450 kyr BP.

Quantitative estimates of the integrated heat content of the entire ocean are required to answer these questions. However, obtaining a representative estimate of the whole-ocean heat content represents a formidable challenge (Hoffman et al., 2017; Shakun et al., 2015). Marine sediments are available for many but not all ocean areas and record only local sea surface conditions or deep-ocean temperatures at the coring site. Moreover, the proxy temperature information gained from marine sediments may be affected by biological processes, sea level changes, or limited precision, especially for cold deep-ocean temperatures. Finally, obtaining globally representative MOT from individual marine archives requires the compilation of a globally representative set of marine sediment cores, including rigorous cross-dating of all the records, which is difficult for the pre-¹⁴C period.

The measurement of atmospheric noble gases trapped in Antarctic ice provides a unique opportunity to reconstruct changes in MOT independently from marine records (Craig and Wiens, 1996; Bereiter et al., 2018b; Headly and Severinghaus, 2007; Ritz et al., 2011) and may overcome these limitations. Changes in the atmospheric concentrations of xenon, krypton, and nitrogen (we liberally add N₂ to the noble gases, as its atmospheric residence time amounts to tens of millions of years) mirror the anomalies in total ocean heat content, as they are driven by changes in their temperature-dependent physical solubilities in ocean water. This allows us to reconstruct MOT in the past from the noble gas composition of a single ice core sample. A minor exception is the small amount of geothermal heating of ocean water, which is not connected to a noble gas change. Note that the MOT at a given time in the past represents a snapshot of the heat (and noble gas) content of all water parcels of the ocean and thus integrates over water bodies with different ventilation ages. Importantly, MOT is therefore not directly linked to mean global sea surface temperature at the same time, a point that has caused some confusion in the past.

Here we extend previous efforts in reconstructing MOT using noble gases in ice cores (Bereiter et al., 2018b; Baggenstos et al., 2019; Shackleton et al., 2020), which were limited to the last and penultimate glacial terminations, by reconstructing the glacial–interglacial MOT range for snapshots of all glacial and interglacial intervals over the last 700 000 years using the EPICA Dome C (EDC) ice core. To obtain quantitative MOT values, large corrections of the raw data for physical transport processes in the firn have to be applied, which affect the precision and accuracy of the MOT data. Moreover, gas fractionation processes in the ice between bubbles and clathrates in connection with potential post-coring gas loss have been observed (Bender et al., 1995), which have to be unambiguously identified to allow us to derive an unbiased MOT estimate, as described in detail in this study. Despite these corrections, we are able to derive quantitative MOT values from high-precision measurements of noble gas mixing ratios in the EDC ice core.

The paper is organized as follows: in Sect. 2 we describe the overall analytical procedure to obtain noble gas isotopic and elemental ratios with the associated uncertainty as well as the corrections that have to be applied to correct for systematic transport effects in the firn column. As EDC is a very low-accumulation site, which exhibits a permanent firn temperature gradient larger than 1 ^∘C, thermal diffusion has to be precisely quantified. Moreover, we describe in detail for the first time the correction of kinetic fractionation by non-diffusive transport, adding an additional layer of complexity and uncertainty. Another process acting on the noble gas composition in the EDC ice core is post-coring gas loss that affects samples from the bubble–clathrate transition zone and, more severely, fully clathrated ice from the deepest ice at Dome C, which is close to the pressure melting point. These effects are described in detail for the first time in Sect. 3, and a hypothesis on how these lead to systematic noble gas fractionation is presented. For the reader only interested in the final results of our MOT reconstruction over the last 700 kyr and their discussion, we refer directly to Sect. 4 and the Conclusions section.

2 Measurement and corrections

2.1 Determination of noble gas elemental and isotopic ratios in ice cores

The data set presented here consists of 88 EDC ice samples (39 samples span the last 40 kyr, and 49 samples are from MIS 5 to MIS 17; Baggenstos et al., 2019) that were analyzed at the University of Bern. In order to test the representativeness of our EDC MOT data, 4 samples from the EPICA Dronning Maud Land (EDML) ice core, 4 samples from the Greenland Ice Core Project (GRIP) ice core, 4 samples from the Greenland NEEM ice core, 4 samples from the Antarctic Talos Dome ice core, and 13 blue ice samples from Taylor Glacier were measured. The Taylor Glacier samples stem from a lab intercomparison with the Scripps Institution of Oceanography in San Diego and showed excellent agreement within the analytical uncertainties between the two labs. Two samples from EDC (age > 700 kyr BP) were discarded due to bad ice quality together with drill fluid contamination, and one sample was rejected due to a procedural problem.

The gas extraction and processing broadly follow the method described in detail by Bereiter et al. (2018a). In short, after shaving off the outer surface ∼600 g of ice is melt-extracted (equivalent to ∼55 mL at the standard temperature and pressure (STP) of extracted air) and the released air is quantitatively frozen into stainless-steel dip tubes in a helium exchange cryostat. To equilibrate the sample and to avoid fractionation by thermal gradients in the dip tube, the dip tubes are kept for at least 10 h in a ventilated isothermal box. The samples are then split into two aliquots. The larger (∼53 mL STP) aliquot is exposed to a $\mathrm{Zr}/\mathrm{Al}$ getter alloy at 900 ^∘C to remove all reactive gases and is analyzed in a peak-jumping routine on a Thermo Finnigan MAT 253 dual-inlet isotope ratio mass spectrometer (IRMS) for the isotopic ratios of xenon, krypton, and argon as well as their elemental ratios. The smaller aliquot (∼2 mL STP) is passed through a CO₂ trap and measured on a Thermo Finnigan MAT Delta V dual-inlet IRMS in parallel for argon, oxygen, and nitrogen isotopes and their elemental ratios. The measurements were corrected for pressure imbalance and chemical slope according to the procedure described in Severinghaus et al. (2003). All data are reported in the delta notation with respect to a modern atmosphere standard collected outside the lab in Bern, Switzerland. Note that this atmospheric value reflects the current dissolution equilibrium of the studied gases between the atmosphere and the global ocean; this atmospheric standard is not subject to fractionation processes in the firn that affect air samples from ice cores. Accordingly, in the absence of any long-term changes in the total abundance of noble gas isotopes in the combined atmosphere and ocean reservoirs, each determined ice core value provides a measure of the difference in the ocean heat content at the age of the sample relative to today.

Table 1 Long-term reproducibility of the analyses (1 standard deviation) derived from 51 measurements of outside air on the MAT 253 and 46 measurements of outside air on the Delta V Plus.

The long-term reproducibility of our system, which was determined using outside air samples, may be found in Table 1. When normalized to one unit of mass difference, all isotopic ratios can be quantified to better than 10 per meg with uncertainties as low as 5 per meg per difference mass unit for ${\mathit{\delta }}^{\mathrm{136}/\mathrm{129}}\mathrm{Xe}$, 7 per meg per difference mass unit for ${\mathit{\delta }}^{\mathrm{86}/\mathrm{84}}\mathrm{Kr}$, and as low as 1.5 per meg per unit mass difference for ${\mathit{\delta }}^{\mathrm{40}/\mathrm{36}}\mathrm{Ar}$. Only δ¹⁵N values measured in this study did not yet reach the expected precision per unit mass difference. Accordingly, we refrain from using δ¹⁵N to correct for gravitational or thermal diffusion fractionation in this study.

Despite this high analytical long-term reproducibility derived from outside air samples, some species showed a much higher scatter in ice core samples due to remnants of drill fluid in the gas sample after processing. Drill fluid interference in the mass spectrometer alters δ¹⁵N and ${\mathit{\delta }}^{\mathrm{40}/\mathrm{36}}\mathrm{Ar}$ in the analysis of the smaller, non-gettered aliquot. The affected samples can be identified using an outlier routine and removed from the data set. Nitrogen isotopes were precisely measured by Dreyfus et al. (2010), Landais et al. (2013), and Bréant et al. (2019) in previous studies. Despite the lower precision of our δ¹⁵N analyses, our results for samples not affected by drill fluid contamination are on average in good agreement with these values (see also Fig. 3). Nevertheless, we refrain from using the δ¹⁵N and ${\mathit{\delta }}^{\mathrm{40}/\mathrm{36}}\mathrm{Ar}$ values measured using the smaller aliquot.

Unexpectedly, we also see outliers in ${\mathit{\delta }}^{\mathrm{86}/\mathrm{82}}$Kr and ${\mathit{\delta }}^{\mathrm{84}/\mathrm{82}}$Kr values in the gettered larger aliquot in the same samples in which we find the δ¹⁵N and ${\mathit{\delta }}^{\mathrm{40}/\mathrm{36}}\mathrm{Ar}$ outliers, which points to an interference at mass 82. We have no conclusive evidence yet to indicate what is causing this interference. Potentially, the large variety of higher organic compounds in the drill fluid allows some component not to be completely gettered if the drill fluid contamination is too large and causes the mass 82 interference. Further lab experiments are required to elucidate this issue. In summary, we only use ${\mathit{\delta }}^{\mathrm{86}/\mathrm{84}}\mathrm{Kr}$ values in our further data evaluation.

To derive MOT, we need to know the changes in the noble gas mixing ratios ($\mathit{\delta }\mathrm{Kr}/{\mathrm{N}}_{\mathrm{2}}$, $\mathit{\delta }\mathrm{Xe}/{\mathrm{N}}_{\mathrm{2}}$, and $\mathit{\delta }\mathrm{Xe}/\mathrm{Kr}$). To obtain them, we mathematically combine $\mathit{\delta }\mathrm{Xe}/\mathrm{Ar}$, $\mathit{\delta }\mathrm{Kr}/\mathrm{Ar}$, and $\mathit{\delta }\mathrm{Ar}/{\mathrm{N}}_{\mathrm{2}}$. Ar itself is not used as a reference element because argon is preferentially excluded relative to N₂, xenon, and krypton during bubble formation at the firn–ice transition (Severinghaus and Battle, 2006). The uncertainties calculated by error propagation of the analytical errors in the elemental ratios used in this calculation are 0.136 ‰ and 0.181 ‰ for $\mathit{\delta }\mathrm{Kr}/{\mathrm{N}}_{\mathrm{2}}$ and $\mathit{\delta }\mathrm{Xe}/{\mathrm{N}}_{\mathrm{2}}$, respectively.

Figure 1 Measured delta values of the EDC elemental ratios from 700 kyr BP to the present. The analytical uncertainties of the elemental ratios are too small (on the order of 0.1‰; see text) to be displayed here. The red shaded intervals indicate interglacials as identified by Masson-Delmotte et al. (2010).

2.2 Inferring atmospheric noble gas ratios from raw data

While the troposphere is well-mixed through turbulent processes, the low permeability of the firn restricts bulk flow; thus, the gas transport is controlled by molecular diffusion. All measured isotope and elemental ratios we extract from the ice core samples (reported as “raw ratios”) are therefore highly fractionated with respect to the atmospheric values (as illustrated by the strongly positive δ values in Fig. 1), mainly due to gravitational enrichment in the diffusive zone of the firn but also thermal diffusion if a temperature gradient is present in the firn column (Schwander et al., 1988) or non-diffusive transport processes that disturb the diffusive equilibrium (so-called kinetic fractionation). These systematic fractionations in the isotopic and elemental ratios have to be corrected for, as illustrated in Fig. 2.

Figure 2 Flow scheme to infer paleo-atmospheric elemental ratios and finally MOT using the measured elemental ($\mathit{\delta }\mathrm{Kr}/{\mathrm{N}}_{\mathrm{2}}$, $\mathit{\delta }\mathrm{Xe}/{\mathrm{N}}_{\mathrm{2}}$, and $\mathit{\delta }\mathrm{Xe}/\mathrm{Kr}$) and isotopic ratios (${\mathit{\delta }}^{\mathrm{40}/\mathrm{36}}\mathrm{Ar}$, ${\mathit{\delta }}^{\mathrm{86}/\mathrm{84}}\mathrm{Kr}$, ${\mathit{\delta }}^{\mathrm{136}/\mathrm{129}}\mathrm{Xe}$). Values are indicated by dark grey square boxes, while applied corrections are illustrated by rounded boxes in lighter grey. In a first step, the ${\mathit{\delta }}^{\mathrm{40}/\mathrm{36}}\mathrm{Ar}$ isotope ratios have to be corrected for changes in ⁴⁰Ar due to geological outgassing. The isotopic and elemental ratios are then corrected for thermal diffusion and kinetic fractionation. The temperature difference between the top and the bottom of the diffusive column height required for the correction of thermal diffusion is derived using either a firn column model (Michel, 2016) or the full matrix of noble gas isotopes. For the kinetic fractionation, we assume fixed ratios between isotopic and elemental ratios (Birner et al., 2018). Finally, the corrected isotopic ratios are used to correct the elemental ratios for the gravitational enrichment. The corrected elemental ratios are then translated into MOT using an ocean atmosphere box model. Here, we correct for the influences of sea level (SL) change and the saturation state (SAT) of the global ocean. In a last step, the resulting MOTs are corrected for their Holocene offset (see text for details). The approximate range of the respective correction performed on all our samples for MOT values can be found in the colored boxes. Positive values (red boxes) imply that the correction leads to an increase in absolute MOT values. It is important to stress that the size of the correction is not a measure of the uncertainty associated with it. The gravitational enrichment correction, for example, dominates all other corrections by orders of magnitude. Nonetheless, the analytical uncertainty of that correction is only minor, as it is largely given by the uncertainty in diffusive column height, the exact determination of which is one of our analytical foci. Also note that it is the size of the range of a correction, rather than the absolute value of the correction, that has an influence on glacial–interglacial MOT differences. For example, although the correction of elemental ratios for kinetic fractionation appears to be the second most influential correction, when considering its absolute value, its range is only narrow. Thus, it alters the glacial–interglacial MOT difference only slightly.

Moreover, post-coring gas losses accompanied by gas fractionation processes have been observed in our samples in the bubble–clathrate transition zone (BCTZ) and in poorly preserved and cracked ice, mainly in the deepest section of the ice core. While the gas fractionation processes mentioned above can be corrected for using the full matrix of noble gases and their isotopes, the gas loss fractionation cannot yet be quantitatively corrected for, and affected samples have to be unambiguously identified using clear detection criteria and removed from the data set. In the following sections, the systematic correction of the isotopic and elemental ratios is described, while the post-coring gas loss and its detection are described in Sect. 3.

2.2.1 Systematic processes altering the isotopic and elemental ratios

Geological outgassing of ⁴⁰Ar

In contrast to the nitrogen, krypton, and xenon isotopes as well as ³⁶Ar and ³⁸Ar, the atmospheric abundance of ⁴⁰Ar has gradually increased over time. This atmospheric ⁴⁰Ar increase is caused by integrated ⁴⁰Ar outgassing throughout Earth's history due to the radioactive decay of ⁴⁰K in the crust and mantle (Bender et al., 2008). Assuming the ⁴⁰Ar increase is constant in time, the effect can be corrected if the age of the sample is known. In this study, we use ${\mathit{\delta }}^{\mathrm{40}/\mathrm{36}}\mathrm{Ar}$ values corrected for outgassing according to Bender et al. (2008).

$$\begin{array}{}\text{(1)}& \mathrm{\Delta }{\mathit{\delta }}^{\mathrm{40}/\mathrm{36}}{\text{Ar}}_{\text{outgas}}=\mathrm{6.6}\times {\mathrm{10}}^{-\mathrm{5}}\mathrm{‰}/\text{kyr}\times t_{\text{gas}}\end{array}$$

Here, the gas age t_gas has units of thousands of years before present (kyr BP). The outgassing effect on MOT is small for samples from the last transition (<0.1 ^∘C), but it becomes more important for samples in the deeper and older ice (on the order of 1 ^∘C for the oldest samples).

Gravitational enrichment

Gas transport in the pore space of the firn column is controlled by molecular diffusion below a convective zone in the top meters of the firn. In this diffusive column molecules of different masses are separated, leading to gravitational enrichment of the heavier gases and isotopes at the bottom of the firn column where bubble close-off occurs. Thus, all elemental and isotopic ratios are strongly enriched in the heavier gases compared to the atmosphere according to

$$\begin{array}{}\text{(2)}& \begin{array}{rl}& \mathrm{\Delta }{\mathit{\delta }}_{\text{grav}}\left(\text{z}\right)={\mathit{\delta }}_{\text{grav}}\left(\text{z}\right)-{\mathit{\delta }}_{\text{atm}}=\left(\right({\mathit{\delta }}_{\text{atm}}+\mathrm{1})\\ & \times e^{\frac{\mathrm{\Delta }mgz}{RT}}-\mathrm{1})-{\mathit{\delta }}_{\text{atm}},\end{array}\end{array}$$

where g is the local gravitational acceleration, R the gas constant, and T the mean firn temperature (Craig et al., 1988; Schwander, 1996). To correct for the gravitational enrichment of the elemental ratios at the lock-in depth, we have to know the diffusive column height (DCH) z of the firn (where z is defined positive downward), which can be directly derived from the isotopic enrichment of the individual gas species after corrections for thermal diffusion and kinetic effects (see below).

Thermal diffusion

The second diffusive process leading to fractionation of the air in the firn column is thermal diffusion. It refers to the fractionation of a gas mixture in the presence of a temperature gradient. This leads to an enrichment of the heavier components at the cold end of the firn column according to the laboratory-determined thermal diffusion sensitivity Ω^X/Y for two gas or isotope species X and Y, where Ω^X/Y is specific for the gas mixture (Severinghaus et al., 2003).

$$\begin{array}{}\text{(3)}& \mathrm{\Delta }{\mathit{\delta }}_{\text{therm}}={\mathrm{\Omega }}^{\text{X/Y}}\times \mathrm{\Delta }T\end{array}$$

Here, ΔT is defined as the temperature difference between the top and the bottom of the DCH. As EPICA Dome C is a low-accumulation site, the mean annual firn temperature increases with depth (in the absence of temporal climate changes at the surface; Ritz et al., 1982) due to the geothermal heat flux at the bottom, and the gas enclosed in bubbles is expected to be slightly depleted by thermal diffusion relative to its gravitational value. Note that this thermal fractionation is gas-species-dependent and much larger for the lighter gases. Here we use the thermal diffusion sensitivities given by Kawamura et al. (2013) and Headly (2008).

Kinetic fractionation: the heavy isotope deficit

If only diffusive processes occurred in the firn column, the differences between any of the measured isotope ratios could be used to unambiguously reconstruct the thermal and gravitational fractionation components using the well-known thermal diffusion sensitivity parameters (Kawamura et al., 2013; Headly, 2008). In particular, as the gravitational enrichment of two gases only differs according to the mass difference between two gases, the gravitational enrichment normalized to unit mass difference should be the same for all gases. However, even after the thermal diffusion correction the isotopic ratios of different gases per unit mass difference still reveal systematic offsets (see Fig. 3a), indicating that there is another mechanism at play that we have yet to correct for.

The reason for the systematic differences in the isotopic ratios per unit mass difference between the different gases at Dome C in Fig. 3 is the occurrence of advective processes in the firn such as turbulent mixing at the surface (Kawamura et al., 2013), barometric pumping, and the net vertical movement of air by ice flow and compression (Buizert and Severinghaus, 2016; Birner et al., 2018) that all disturb the diffusive equilibrium in the firn column. The term coined for this difference in measured isotopic ratios per unit mass difference after thermal correction is heavy isotope deficit (HID) (Buizert and Severinghaus, 2016) or differential kinetic fractionation (Birner et al., 2018) between the isotopic or elemental ratio X and another isotopic or elemental ratio Y. Further on, we will use the difference between the kinetic fractionation of an isotopic or elemental ratio X per unit mass difference relative to the kinetic fractionation in N₂ isotopes (Kawamura et al., 2013; Birner et al., 2018; Buizert and Severinghaus, 2016). For example, for Ar isotopes this fractionation is defined as

$$\begin{array}{}\text{(4)}& {\mathit{ϵ}}_{\text{Ar}-{\text{N}}_{\mathrm{2}}}^{\ast }=\mathrm{\Delta }{\mathit{\delta }}^{\mathrm{40}/\mathrm{36}}{\text{Ar}}_{\text{kin}}^{\ast }-\mathrm{\Delta }{\mathit{\delta }}^{\mathrm{15}}{\text{N}}_{\text{kin}}^{\ast }.\end{array}$$

Here, ^∗ indicates that the values are normalized to unit mass difference, for example ${\mathit{\delta }}^{\mathrm{40}/\mathrm{36}}{\text{Ar}}_{\text{kin}}^{\ast }={\mathit{\delta }}^{\mathrm{40}/\mathrm{36}}{\text{Ar}}_{\text{kin}}/\mathrm{4}$. Note that the heavier (and slowly diffusing) gases like xenon and krypton will be further from diffusive equilibrium than lighter (and faster-diffusing) gases like argon and nitrogen. Another way to put this is that the effective DCH of different gases is not the same, with z being a little shorter for heavier, less diffusive gases.

In our measurements the total deficit of xenon and krypton isotopes relative to argon isotopes per unit mass difference over the last 40 kyr is on average −0.046 ± 0.007 ‰ and −0.029 ± 0.007 ‰, respectively (here and throughout the paper all uncertainties refer to 1σ). This total deficit comes about through the combined effect of different thermal diffusion sensitivities of the different gases and gas-specific kinetic fractionation. Separation of these two effects on the total deficit requires corrections for the thermal gradient, as described in Sect. 2.2.2.

Birner et al. (2018) discovered that while the absolute magnitude of kinetic fractionation $\mathrm{\Delta }{\mathit{\delta }}_{\text{kin}}^{\ast }$ may vary substantially between different ice core sites, the kinetic fractionations between two gases occur in fixed ratios between isotopic and elemental ratios for a variety of firn regimes. Thus, although we do not know the absolute kinetic fractionation at Dome C for an individual gas for past firn conditions a priori, we can use the results by Birner et al. (2018) together with our multiple gas species and transient firn temperature modeling to correct for the HID and other gas transport effects. It is worth noting that the ratios of the mean kinetic fractionations over the last 40 kyr as displayed in Fig. 3 agree well within uncertainties with the ratios predicted by Birner et al. (2018). Note that the kinetic fractionations for the last 40 kyr displayed in Fig. 3 also suggest, apart from substantial analytical scatter, some systematic variations. As these systematic variations are generally in line with the fractionation ratios predicted by Birner et al. (2018), the changes are likely caused by changing kinetic fractionation conditions with time (stronger barometric pumping, higher wind speeds at the surface, etc.). However, some samples older than 40 kyr show much higher deviations from the mean that cannot be explained by analytical uncertainty or changing kinetic fractionation. Accordingly, we do not use the individual HID values of each sample as displayed in Fig. 3b for samples older than 40 kyr to quantify ${\mathit{ϵ}}_{\text{Ar}-{\text{N}}_{\mathrm{2}}}^{\ast }$. Instead, we use the mean isotopic fractionation derived for the last 40 kyr and assume that this mean value is also representative for the correction of our samples older than 40 kyr. As the choice of ${\mathit{ϵ}}_{\text{Ar}-{\text{N}}_{\mathrm{2}}}^{\ast }$ has a systematic effect on the final reconstructed MOT, we also take the variation of ${\mathit{ϵ}}_{\text{Ar}-{\text{N}}_{\mathrm{2}}}^{\ast }$ over the last 40 kyr into account to quantify the systematic uncertainty introduced by this choice. Thus, we use the mean kinetic fractionation over the last 40 kyr plus (minus) its standard deviation as upper (lower) bounds of the kinetic fraction to calculate MOT. We regard this as a systematic uncertainty in our MOT reconstruction and separate this uncertainty from the stochastic uncertainty introduced by the analytical error.

Gas loss effects on isotopic ratios

As mentioned in the previous paragraph, ice core samples with heavy post-coring gas loss exhibit enriched ${\mathit{\delta }}^{\mathrm{40}/\mathrm{36}}\mathrm{Ar}$, while δ¹⁵N is unaffected (J. Severinghaus, personal communication, 2021). Due to their large size, the krypton and xenon isotopes are also not affected. Correction coefficients for argon isotope enrichment in bubble ice have been suggested by Kobashi et al. (2008b), but it is not certain whether these correction factors describe the gas loss quantitatively and in particular whether they can be used for ice from any depth and for all ice cores. We refrain from trying to correct this effect because we do not see evidence for an isotope fractionating effect by artifactual gas loss (see Sect. 3.1) and because no correction coefficients have been established for clathrate ice.

Figure 3 (a) Isotopic enrichment of the EDC samples from 700 kyr BP to the present. All isotopic data are shown per unit mass difference and reported relative to modern air sampled in Bern, Switzerland. The data are corrected for thermal fractionation as described in the text. In the case of δ¹⁵N, data points contaminated with drilling fluid have been excluded. The circles in darker green, olive, and lighter green represent measurements from the EDC ice core from Dreyfus et al. (2010), Landais et al. (2013), and Bréant et al. (2019), which have been corrected for thermal fractionation using our model-derived firn temperatures. (b) Difference between the isotopic ratios of the EDC samples after correction for thermal fractionation. The dashed line is the mean difference from the samples over the last 40 kyr, with the shaded areas representing their standard deviation. The red shaded intervals and the MIS numbers on top indicate interglacials as identified by Masson-Delmotte et al. (2010).

2.2.2 Correcting isotopic ratios for transport processes in the firn column

As outlined above, the measured isotopic ratio ${\mathit{\delta }}_{\text{meas}}^{\ast }$ is influenced by a series of processes that alter the paleo-atmospheric value ${\mathit{\delta }}_{\text{atm}}^{\ast }$. We can write the following for the measured isotopic ratio.

$$\begin{array}{}\text{(5)}& {\mathit{\delta }}_{\text{meas}}^{\ast }={\mathit{\delta }}_{\text{atm}}^{\ast }+\mathrm{\Delta }{\mathit{\delta }}_{\text{grav}}^{\ast }+\mathrm{\Delta }{\mathit{\delta }}_{\text{therm}}^{\ast }+\mathrm{\Delta }{\mathit{\delta }}_{\text{kin}}^{\ast }\end{array}$$

Note that we assume ocean solubility effects on the isotopic ratios to be negligible (in contrast to the elemental ratios, for which we use these effects to quantify MOT). This can be justified using a two-box model calculation with the recently determined solubility fractionation factors in seawater (Seltzer et al., 2019) and imposed ocean temperature differences, which yields glacial–interglacial changes in atmospheric noble gas isotope composition of less than 0.001 ‰ per unit mass difference. We can thus safely assume that for isotopic ratios ${\mathit{\delta }}_{\text{atm}}^{\ast }=\mathrm{0}$ (where for Ar isotopes the atmospheric value has to be corrected for Ar outgassing as described in Sect. 2.2.2).

Using the measured isotopic ratios we are able to quantify the fractionation effects and correct the isotopic and elemental ratios. Here we explore two different correction pathways that differ in the derivation of ΔT, which is required to quantify the thermal diffusion effect. Both possibilities use the matrix of isotopic ratios to constrain the gravitational enrichment and the HID. However, one possibility uses a 1D ice flow model connected to a dynamic heat advection and diffusion model (Michel, 2016) to quantify the thermal diffusion effect, while the other uses the full matrix of noble gas isotopes measured in this study to quantify ΔT empirically.

(a) Correction using a firn temperature model

To correct for the strong gravitational enrichment (large mass difference) we can use the precisely measured isotopic ratios after correcting for thermal diffusion and kinetic fractionation. Due to the firn temperature slightly increasing with depth at EDC and due to the negative kinetic fractionation for heavier gases, the corrected isotopic ratios are slightly increased compared to the measured ratio. The results for the thermally corrected isotopic ratios reflecting the remnant kinetic fractionation are shown in Fig. 3a, where ΔT used for the thermal correction of each individual sample was calculated using an ice flow–heat flow model as described below. The isotopic ratio least affected by kinetic fractionation is δ¹⁵N; however, our δ¹⁵N measurements are not as precise as ${\mathit{\delta }}^{\mathrm{40}/\mathrm{36}}\mathrm{Ar}$ per unit mass difference, and some of the δ¹⁵N values had to be discarded due to drill fluid contamination affecting the mass spectrometric analysis. Accordingly, we use ${\mathit{\delta }}^{\mathrm{40}/\mathrm{36}}\mathrm{Ar}$ as our reference isotope. In many applications to reconstruct temperature using isotopic ratios of permanent gases (e.g., Kobashi et al., 2008a), the small kinetic fraction of Ar relative to N₂ is neglected and the gravitational enrichment per unit mass difference $\mathrm{\Delta }{\mathit{\delta }}_{\text{grav}}^{\ast }$ in units of per mill is calculated according to

$$\begin{array}{}\text{(6)}& \mathrm{\Delta }{\mathit{\delta }}_{\text{grav}}^{\ast }={\mathit{\delta }}^{\mathrm{40}/\mathrm{36}}{\text{Ar}}^{\ast }-{\mathrm{\Omega }}^{\mathrm{40}/\mathrm{36}\ast }\times \mathrm{\Delta }T,\end{array}$$

where the ^∗ in ${\mathrm{\Omega }}^{\mathrm{40}/\mathrm{36}\ast }$ indicates that the thermal sensitivity has been normalized to unit mass difference.

However, firn gas pumping and modeling experiments (Kawamura et al., 2013) show that all gases are subject to kinetic fractionation, and therefore even Ar isotopes show a small heavy isotope deficit ${\mathit{ϵ}}_{\text{Ar}-{\text{N}}_{\mathrm{2}}}^{\ast }$ after correction of thermal diffusion effects (Buizert and Severinghaus, 2016) compared to nitrogen isotopes. Thus, $\mathrm{\Delta }{\mathit{\delta }}_{\text{grav}}^{\ast }$ should be calculated according to

$$\begin{array}{}\text{(7)}& \begin{array}{rl}& \mathrm{\Delta }{\mathit{\delta }}_{\text{grav}}^{\ast }={\mathit{\delta }}^{\mathrm{40}/\mathrm{36}}{\text{Ar}}^{\ast }-{\mathrm{\Omega }}^{\mathrm{40}/\mathrm{36}\ast }\cdot \mathrm{\Delta }T-\mathrm{\Delta }{\mathit{\delta }}^{\mathrm{40}/\mathrm{36}}{\text{Ar}}_{\text{kin}}^{\ast }\\ & ={\mathit{\delta }}^{\mathrm{40}/\mathrm{36}}{\text{Ar}}^{\ast }-{\mathrm{\Omega }}^{\mathrm{40}/\mathrm{36}\ast }\cdot \mathrm{\Delta }T-{\mathit{ϵ}}_{\text{Ar}-{\text{N}}_{\mathrm{2}}}^{\ast }-\mathrm{\Delta }{\mathit{\delta }}^{\mathrm{29}/\mathrm{28}}{\text{N}}_{\mathrm{2},\text{kin}}^{\ast }.\end{array}\end{array}$$

The absolute kinetic fractionation $\mathrm{\Delta }{\mathit{\delta }}^{\mathrm{40}/\mathrm{36}}{\text{Ar}}_{\text{kin}}^{\ast }$ in Eq. (7) is not known a priori and is dependent on firn and meteorological conditions at the site that may have changed in the past. Using the equivalent equation for Xe isotopes and the linear relationship of the kinetic fractionations given by Birner et al. (2018), ${\mathit{ϵ}}_{\text{Ar}-{\text{N}}_{\mathrm{2}}}^{\ast }$ can be calculated from our measurements according to

$$\begin{array}{}\text{(8)}& \begin{array}{rl}& {\mathit{ϵ}}_{\text{Ar}-{\text{N}}_{\mathrm{2}}}^{\ast }=({\mathit{\delta }}^{\mathrm{136}/\mathrm{129}}{\text{Xe}}_{\text{meas}}^{\ast }-{\mathit{\delta }}^{\mathrm{40}/\mathrm{36}}{\text{Ar}}_{\text{meas}}^{\ast }\\ & -({\mathrm{\Omega }}^{\mathrm{136}/\mathrm{129}\ast }-{\mathrm{\Omega }}^{\mathrm{40}/\mathrm{36}\ast })\cdot \mathrm{\Delta }T)/\mathrm{6.3}.\end{array}\end{array}$$

Note that while including the kinetic fractionation may improve the accuracy of the reconstruction, it slightly decreases the precision, as the analytical error for Xe isotopes also comes into play.

ΔT in Eqs. (7) and (8) is unknown but can be estimated for each individual sample using an ice flow–heat flow model that solves the ice advection–firn compaction and the heat advection and diffusion equation transiently over the last million years. Here we used the model by Michel (2016) and a parameterization of heat diffusion and heat conductivity according to Cuffey and Paterson (2010). The density of the firn was estimated by the Herron–Langway model (Michel, 2016). Heat diffusivity of the very porous top 10 m of the firn was corrected according to Weller and Schwerdtfeger (1970). The recent firn temperature difference ΔT in this run agreed within ±0.2 ^∘C between the model and measured firn temperature profiles at Dome C (Buizert et al., 2020). We assume that ±0.2^∘C represents the analytical uncertainty of the model. Surface temperature and snow accumulation used in the model are prescribed using the EDC temperature record over the last 800 000 years (Jouzel et al., 2007; Bazin et al., 2013), and a scaled version of the stacked marine δ¹⁸O record (Lisiecki and Raymo, 2005) is used for earlier times. For an optimal fit with the measured temperature profile over the entire ice thickness at Dome C, a +2 ^∘C correction of glacial temperatures was applied by Michel (2016); however, this has only a very small effect on the firn temperature difference. Ice thickness changes over the last 800 kyr are taken from Parrenin et al. (2007), and a scaled version of the stacked marine δ¹⁸O record is used for earlier times. The constant geothermal heat flux and the thinning function were chosen to optimize the agreement of the modeled temperature profile with borehole temperature measurements (Buizert et al., 2020) and the age profile at Dome C (Bazin et al., 2013).

The DCH affects both the gravitational enrichment and the absolute value of the temperature gradient. Accordingly, we use an iterative approach for Eq. (7) to determine the firn temperature gradient from the model results. In the first step we neglect the temperature gradient, and thus thermodiffusion effects, to estimate the DCH:

$$\begin{array}{}\text{(9)}& z=\text{ln}({\displaystyle \frac{\mathrm{\Delta }{\mathit{\delta }}_{\text{grav}}^{\ast }}{\mathrm{1000}}}+\mathrm{1})\cdot {\displaystyle \frac{\text{R}T}{\mathrm{\Delta }{m}^{\ast }\text{g}}},\end{array}$$

where Δm^∗ is the unit mass difference in kilograms per mole (kg mol⁻¹). After this first step we read out the modeled firn temperature difference for this depth and use it for the thermal diffusion correction in the next step. Using this iterative approach the firn column depth converges after two iterations. This model-based analysis shows that the firn column depth is about 70–80 m for glacial periods and 90–100 m for interglacials.

The kinetic fractionation ${\mathit{ϵ}}_{\text{Ar}-{\text{N}}_{\mathrm{2}}}^{\ast }$ calculated using modeled ΔT is on average $-\mathrm{0.009}\pm \mathrm{0.001}$ ‰ (1σ) over the last 40 kyr (note that the HID values given in Fig. 3b represent a measure of ${\mathit{ϵ}}_{\text{Ar}-{\text{N}}_{\mathrm{2}}}^{\ast }$ scaled by a factor of 6.3, 4.25, and 2.05 for ${\mathit{ϵ}}_{\text{Xe}-\text{Ar}}^{\ast }$, ${\mathit{ϵ}}_{\text{Kr}-\text{Ar}}^{\ast }$, and ${\mathit{ϵ}}_{\text{Xe}-\text{Kr}}^{\ast }$, respectively, as derived from Birner et al., 2018). As described above, we use the mean plus or minus the standard deviation of the kinetic fractionation over the last 40 kyr to also correct the samples older than 40 kyr and to assess the uncertainty introduced by this correction.

The last unknown in Eqs. (7) and (4) is the amount of kinetic fractionation for nitrogen isotopes $\mathrm{\Delta }{\mathit{\delta }}^{\mathrm{15}}{\text{N}}_{\text{kin}}^{\ast }$. Modeling and firn air studies by Buizert and Severinghaus (2016) and Birner et al. (2018) show that at WAIS Divide, i.e., a site on the West Antarctic Plateau, this fractionation is less than 0.005 ‰, while at Law Dome, a site subject to strong barometric pumping, it may be larger (J. Severinghaus, personal communication, 2021). Assuming a kinetic fractionation at Dome C that is similar to WAIS, this small fractionation in $\mathrm{\Delta }{\mathit{\delta }}^{\mathrm{15}}{\text{N}}_{\text{kin}}^{\ast }$ translates into an offset in the final MOT values on the order of 0.01 ^∘C and can therefore be safely neglected in our MOT reconstruction.

Figure 4 Comparison of the firn temperature difference ΔT_mod, as derived from the ice flow–heat flow model of Michel (2016) (red squares), and ΔT_iso obtained using measured isotopic ratios only (black squares). Dashed lines in red and black indicate the overall mean values for ΔT_mod and ΔT_iso, while the dotted line in cyan represents the mean ΔT_iso when considering samples from the first 40 kyr only. The red shaded intervals and the MIS numbers on top indicate interglacials as identified by Masson-Delmotte et al. (2010).

In the following, we assume $\mathrm{\Delta }{\mathit{\delta }}^{\mathrm{15}}{\text{N}}_{\text{kin}}^{\ast }=\mathrm{0}$ for Dome C, where synoptic pressure variations and wind speeds are rather low. Neglecting potential errors in the thermal diffusion sensitivities, the uncertainty in $\mathrm{\Delta }{\mathit{\delta }}_{\text{grav}}^{\ast }$ introduced by this model-based method to correct the data is 0.003 ‰ (1σ).

(b) Correction using measured isotopic ratios only

Using the firn temperature model as in (a) assumes that the firn temperature difference in the model over the DCH realistically resembles the temperature difference seen by the gas at the age of each individual sample. The model neglects, for example, seasonal variations in the temperature gradient, which may not cancel out in the long-term annual average, and also assumes that the current parameterizations used for firn density, heat conduction, and heat diffusivity profiles in the firn are also representative for past conditions. An alternative way is therefore to use only experimental evidence from our full matrix of isotope ratios. Assuming that gravitational enrichment, thermal diffusion, and kinetic fractionation completely describe the fractionation processes occurring in the firn column, this approach provides a purely measurement-based estimate of ${\mathit{ϵ}}_{\text{Ar}-{\text{N}}_{\mathrm{2}}}^{\ast }$ and ΔT.

The equivalent to Eq. (8) for Kr isotopes is

$$\begin{array}{}\text{(10)}& \begin{array}{rl}& {\mathit{ϵ}}_{\text{Ar}-{\text{N}}_{\mathrm{2}}}^{\ast }=({\mathit{\delta }}^{\mathrm{86}/\mathrm{84}}{\text{Kr}}_{\text{meas}}^{\ast }-{\mathit{\delta }}^{\mathrm{40}/\mathrm{36}}{\text{Ar}}_{\text{meas}}^{\ast }\\ & -({\mathrm{\Omega }}^{\mathrm{86}/\mathrm{84}*}-{\mathrm{\Omega }}^{\mathrm{40}/\mathrm{36}\ast })\cdot \mathrm{\Delta }T)/\mathrm{4.25},\end{array}\end{array}$$

with the constant ratio of kinetic fractionations for different isotopes given by Birner et al. (2018).

Solving Eq. (8) for ΔT and inserting it into Eq. (10), it follows that

$$\begin{array}{}\text{(11)}& \begin{array}{rl}& {\mathit{ϵ}}_{\text{Ar}-{\text{N}}_{\mathrm{2}}}^{\ast }\\ & ={\displaystyle \frac{({\mathit{\delta }}^{\mathrm{86}/\mathrm{84}}{\text{Kr}}^{\ast }-{\mathit{\delta }}^{\mathrm{40}/\mathrm{36}}{\text{Ar}}^{\ast })-\mathit{\alpha }({\mathit{\delta }}^{\mathrm{136}/\mathrm{129}}{\text{Xe}}^{\ast }-{\mathit{\delta }}^{\mathrm{40}/\mathrm{36}}{\text{Ar}}^{\ast })}{\mathrm{4.25}-\mathrm{6.3}\mathit{\alpha }}},\end{array}\end{array}$$

where $\mathit{\alpha }=({\mathrm{\Omega }}^{\mathrm{86}/\mathrm{84}{}^{\ast }}-{\mathrm{\Omega }}^{\mathrm{40}/\mathrm{36}\ast })/({\mathrm{\Omega }}^{\mathrm{136}/\mathrm{129}\ast }-{\mathrm{\Omega }}^{\mathrm{40}/\mathrm{36}\ast })$. Note that the mean kinetic fractionation ${\mathit{ϵ}}_{\text{Ar}-{\text{N}}_{\mathrm{2}}}^{\ast }$ derived using this data-based approach is on average −0.009 ‰ over the last 40 kyr and thus in perfect agreement with the results from the firn temperature model approach.

Using Eq. (11) in Eq. (8) finally leads to

$$\begin{array}{}\text{(12)}& \begin{array}{rl}& \mathrm{\Delta }T\\ & ={\displaystyle \frac{({\mathit{\delta }}^{\mathrm{136}/\mathrm{129}}{\text{Xe}}^{\ast }-{\mathit{\delta }}^{\mathrm{40}/\mathrm{36}}{\text{Ar}}^{\ast })-\frac{\mathrm{6.3}}{\mathrm{4.25}-\mathrm{6.3}\mathit{\alpha }}\left(\right({\mathit{\delta }}^{\mathrm{86}/\mathrm{84}}{\text{Kr}}^{\ast }-{\mathit{\delta }}^{\mathrm{40}/\mathrm{36}}{\text{Ar}}^{\ast })-\mathit{\alpha }({\mathit{\delta }}^{\mathrm{136}/\mathrm{129}}{\text{Xe}}^{\ast }-{\mathit{\delta }}^{\mathrm{40}/\mathrm{36}}{\text{Ar}}^{\ast }\left)\right)}{{\mathrm{\Omega }}^{\mathrm{136}/\mathrm{129}\ast }-{\mathrm{\Omega }}^{\mathrm{40}/\mathrm{36}\ast }}}.\end{array}\end{array}$$

While this method in principle avoids any potential systematic uncertainties in the firn temperature gradient, it substantially increases the uncertainty of the reconstruction, as now the analytical uncertainties of all three noble gas isotopic ratios have to be taken into account for each individual sample. Using Monte Carlo error propagation, the uncertainty in ${\mathit{ϵ}}_{\text{Ar}-{\text{N}}_{\mathrm{2}}}^{\ast }$ derived in this approach is increased to 0.01 ‰ (i.e., roughly 10 times larger than in the model-based approach), and the uncertainty in ΔT becomes about 8 ^∘C. This translates into an error for $\mathrm{\Delta }{\mathit{\delta }}_{\text{grav}}^{\ast }$ of roughly 0.08 ‰ (i.e., almost 30 times higher than in the model-based approach) and finally into an error in the corrected elemental ratios (see Sect. 2.2.3) of several per mill, which is prohibitive for a reliable reconstruction of MOT for individual samples.

The very high scatter in data-derived ΔT is displayed in Fig. 4. The standard deviation of ΔT over the last 40 kyr is 6 ^∘C, which is of the same order as the expected analytical uncertainty in ΔT of 8 ^∘C derived from error propagation. Hence, this large scatter can be explained by analytical noise and shows that any systematic variation in ΔT over the last 40 kyr cannot be reliably quantified using the data-based approach.

Looking at Fig. 4 it should be stressed that despite the large scatter the mean ΔT over the last 40 kyr derived in this data-based approach is in excellent agreement with ΔT derived in the model-based approach. This points to the fact that using our noble gas isotope measurements we are able to accurately (but not sufficiently precisely) reconstruct ${\mathit{ϵ}}_{\text{Ar}-{\text{N}}_{\mathrm{2}}}^{\ast }$ and ΔT over the last 40 kyr and that thermal diffusion and kinetic fractionation fully explain the systematic offsets of the isotopic ratios.

However, for the samples older than 40 kyr the average ΔT is $-\mathrm{7.7}\pm \mathrm{9.7}$ ^∘C, which is physically not in line with the observed or modeled firn temperature differences and also not in line with the results younger than 40 kyr. As we do not expect that any unknown firn fractionation processes were active prior to 40 kyr BP that did not occur over the last 40 kyr, this apparent offset in ΔT either reflects the large scatter of the data or is related to an isotope fractionating Ar loss process during transport, when the samples got relatively warm for a short period of time.

As mentioned above, using the data-based approach to correct the elemental ratios for each individual sample is not a viable solution, as the scatter in reconstructed ΔT, ${\mathit{ϵ}}_{\text{Ar}-{\text{N}}_{\mathrm{2}}}^{\ast }$, and thus MOT is too large. Therefore, we regard the mean ${\mathit{ϵ}}_{\text{Ar}-{\text{N}}_{\mathrm{2}}}^{\ast }$ and ΔT values over the last 40 kyr to be representative for the entire record and use those to calculate $\mathrm{\Delta }{\mathit{\delta }}_{\text{grav}}^{\ast }$. To quantify the systematic uncertainty introduced by this choice of ${\mathit{ϵ}}_{\text{Ar}-{\text{N}}_{\mathrm{2}}}^{\ast }$ and ΔT, we calculate an upper (lower) estimate of ΔT and ${\mathit{ϵ}}_{\text{Ar}-{\text{N}}_{\mathrm{2}}}^{\ast }$ as defined by their mean values over the last 40 kyr plus (minus) their standard deviation. This provides us with a range of values for $\mathrm{\Delta }{\mathit{\delta }}_{\text{grav}}^{\ast }$ for the systematic corrections of the elemental ratios (see Sect. 2.2.3). We refer to this correction as the mean data-based approach further on.

We stress that while such a constant correction with a systematic range of possible values has a significant impact on the absolute reconstructed MOT value, the impact of the correction on the glacial–interglacial difference in MOT is small. Moreover, the glacial–interglacial difference in MOT is very similar whether one uses the model-based approach in (a) or the mean data-based approach in (b). Accordingly, only MOT differences relative to the Holocene are considered, as was also the case in the study by Shackleton et al. (2020).

In summary, the data-based approach shows that we can quantitatively correct the gas-transport-related fractionations that occur in the firn column; however, the precision of the data-based approach is not sufficient for single samples. Accordingly, in the discussion of our final MOT data in Sect. 4 we do not use the data-based approach and rely on the much more precise model-based correction.

2.2.3 Correction of the elemental ratios

Equivalent to isotopic ratios of one gas species, the ratios of two different gas species (here referred to as elemental ratios) have to be corrected for gravitational enrichment. Due to the large mass difference between two gas species, this gravitational correction of elemental ratios is substantial. In addition, thermal diffusion and kinetic fractionation have to be corrected in the elemental ratios. For the $\mathrm{Kr}/{\mathrm{N}}_{\mathrm{2}}$ ratio this can be summarized according to

$$\begin{array}{}\text{(13)}& \begin{array}{rl}& \mathit{\delta }\text{Kr}/{\text{N}}_{\mathrm{2},\text{atm}}=(\mathit{\delta }\text{Kr}/{\text{N}}_{\mathrm{2}\text{,meas}}-{\mathrm{\Omega }}^{\text{Kr}/{\text{N}}_{\mathrm{2}}}\mathrm{\Delta }T\\ & -({\mathit{ϵ}}_{\text{Kr}/{\text{N}}_{\mathrm{2}}-{\text{N}}_{\mathrm{2}}}^{\ast }+\mathrm{\Delta }{\mathit{\delta }}^{\mathrm{15}}{\text{N}}_{\text{kin}}^{\ast })\cdot \mathrm{\Delta }m/\mathrm{\Delta }{m}^{\ast }+\mathrm{1})\\ & \cdot (\mathrm{\Delta }{\mathit{\delta }}_{\text{grav}}^{\ast }+\mathrm{1}{)}^{-\mathrm{\Delta }m/\mathrm{\Delta }{m}^{\ast }}-\mathrm{1},\end{array}\end{array}$$

where again we set $\mathrm{\Delta }{\mathit{\delta }}^{\mathrm{15}}{\text{N}}_{\text{kin}}^{\ast }=\mathrm{0}$, as the effect of $\mathrm{\Delta }{\mathit{\delta }}^{\mathrm{15}}{\text{N}}_{\text{kin}}^{\ast }$ is negligible compared to other uncertainties. The first factor in brackets on the right-hand side represents the measured value of the elemental ratio corrected for thermal diffusion and kinetic fractionation, while the second factor corrects this isotopic ratio for its gravitational enrichment. ΔT, ${\mathit{ϵ}}_{\text{Kr}/{\text{N}}_{\mathrm{2}}-{\text{N}}_{\mathrm{2}}}^{\ast }$, and $\mathrm{\Delta }{\mathit{\delta }}_{\text{grav}}^{\ast }$ can be calculated from the firn temperature model and from the isotopic data as described above using the fixed ratios of kinetic fractionation given by Birner et al. (2018). Note that the exponent $\mathrm{\Delta }m/\mathrm{\Delta }{m}^{\ast }$ is very large (e.g., 56 for ⁸⁴Kr and ²⁸N₂); thus, the gravitational correction is substantial and sensitive to any analytical error in the isotopic ratios as well as in uncertainties in the thermal correction and the kinetic fractionation.

Here we use both the model-based and the mean data-based reconstruction of ΔT and ${\mathit{ϵ}}_{\text{Kr}/{\text{N}}_{\mathrm{2}}-{\text{N}}_{\mathrm{2}}}$ described above and apply these to Eq. (13) and equivalent equations for $\mathrm{Xe}/{\mathrm{N}}_{\mathrm{2}}$ and $\mathrm{Xe}/\mathrm{Kr}$ ratios. The results of these corrections are shown in panel (a) (model-based reconstruction) and (b) (mean data-based reconstruction) of Fig. 5. Note that the whiskers in Fig. 5 represent the stochastic error (1σ) due to error propagation of the analytical uncertainties, while the size of the boxes represents the systematic uncertainties (1σ) introduced by our choice of mean ΔT and ${\mathit{ϵ}}_{\text{Kr}/{\text{N}}_{\mathrm{2}}-{\text{N}}_{\mathrm{2}}}$.

Figure 5 Comparison of the reconstructed atmospheric ratios of $\mathrm{Kr}/{\mathrm{N}}_{\mathrm{2}}$ (orange), $\mathrm{Xe}/{\mathrm{N}}_{\mathrm{2}}$ (purple), and $\mathrm{Xe}/\mathrm{Kr}$ (blue) (a) as derived from the calibrated firn model of Michel (2016) (b) using measured isotopic ratios only. The reconstructed value is drawn as a black line in the middle of each box. The magnitude of the systematic uncertainty is indicated by the size of the colored boxes, and the statistical uncertainties are depicted as whiskers (1σ each). The red shaded intervals and the MIS numbers on top indicate interglacials as identified by Masson-Delmotte et al. (2010).

The first observation is that the three different atmospheric elemental ratios have different glacial–interglacial amplitudes, primarily reflecting the different temperature sensitivities of gas solubility in the ocean as expected. Secondly, the corrected elemental ratios show that although our data are referenced to recent outside air at Bern, reflecting late Holocene ocean temperatures, the Holocene elemental ratios after all corrections are not equal to zero but systematically shifted to positive values.

Holocene elemental ratios larger than zero have also been recently observed by Shackleton et al. (2020) using Taylor Glacier ice, reflecting our insufficient quantitative understanding of fractionation processes in the firn; this requires more dedicated firn gas sampling and firn gas transport modeling studies that realistically resolve seasonal variations in the gas transport, pressure variability, and kinetic fractionation. However, glacial–interglacial differences in our corrected elemental ratios are not sensitive to this overall offset. If we use changes in atmospheric elemental ratios relative to the Holocene mean, both the model-based and the mean data-based approach provide very similar results. For example, using the model-based approach, the corrected last glacial $\mathit{\delta }\mathrm{Xe}/\mathrm{Kr}$ values are on average about 3.0 ‰ lower than the Holocene values, which is in agreement within uncertainties with the results by Bereiter et al. (2018b) using ice from the WAIS Divide ice core.

The final observation pertains to the scatter of the data. Elemental ratios corrected using ΔT and ${\mathit{ϵ}}_{\text{Kr}/{\text{N}}_{\mathrm{2}}-{\text{N}}_{\mathrm{2}}}$ derived using the data-based approach for each individual sample provide no meaningful results, as the analytical scatter introduced from all three isotopic ratios is too large.

Using the model-based approach, we introduce uncertainty through the elemental ratio analysis, the uncertainty in our model-based ΔT, and systematic uncertainty through our use of the mean ${\mathit{ϵ}}_{\text{Kr}/{\text{N}}_{\mathrm{2}}-{\text{N}}_{\mathrm{2}}}$ over the last 40 kyr to correct the kinetic fractionation. This leads to a substantially smaller scatter than in the individual data-based approach. The stochastic errors (1σ) in δKr / N_2,atm, δXe / N_2,atm, and δXe / Kr_atm introduced by the measurement errors are 0.20 ‰, 0.45 ‰, and 0.32 ‰, respectively, in our model-based approach. The systematic uncertainty in this case is 0.23 ‰, 0.70 ‰, and 0.48 ‰, respectively.

Alternatively, we can use the mean data-based correction approach. Using the mean ΔT and ${\mathit{ϵ}}_{\text{Kr}/{\text{N}}_{\mathrm{2}}-{\text{N}}_{\mathrm{2}}}$ over the last 40 kyr leads to the smallest scatter in our atmospheric elemental ratios, as we average over 39 data points for our correction. However, this also implies a large systematic uncertainty that derives from the choice of mean ΔT and ${\mathit{ϵ}}_{\text{Kr}/{\text{N}}_{\mathrm{2}}-{\text{N}}_{\mathrm{2}}}$. Moreover, both parameters may have been subject to systematic changes over time that we cannot decipher using this mean correction. The stochastic errors (1σ) in δKr / N_2,atm, δXe / N_2,atm, and δXe / Kr_atm introduced by the measurement errors are 1.68 ‰, 2.12 ‰, and 0.66 ‰, respectively, in our mean data-based approach. The systematic uncertainty in this case is on average 1.22 ‰, 1.49 ‰, and 0.40 ‰, respectively.

In summary, our different correction pathways lead to similar results in the differences of atmospheric elemental ratios relative to their Holocene value. In the following, we will use only the atmospheric elemental ratios after correction using the model-based approach, which is based on well-understood physical laws of heat conduction–advection and provides elemental ratios with sufficient precision and accuracy to reconstruct glacial–interglacial changes in MOT after correction for gravitational, thermal, and kinetic fractionation.

2.3 Box model to infer the MOT from atmospheric values

To reconstruct MOT values from the corrected elemental ratios, we use the ocean–atmosphere box model described by Bereiter et al. (2018b), except for one slight update, as we use the solubility equations given by Jenkins et al. (2019). The basic assumption in the model is that N₂, Kr, and Xe are conserved in the ocean–atmosphere system and that their distribution between the two reservoirs is controlled by the temperature-dependent dissolution in the ocean water. The model accounts for the effects of changes in ocean salinity, volume, and atmospheric pressure on the oceanic inventories of each of the three gases, which have individual temperature-dependent solubilities. The records of Lambeck et al. (2014) and Spratt and Lisiecki (2016) are used to reconstruct the sea level for the last 27 kyr and the time interval 28–800 kyr BP, respectively. We assume a constant undersaturation of the heavy noble gases in the ocean by 3 % for Xe and 1.5 % for Kr (Hamme and Severinghaus, 2007; Loose et al., 2016; Hamme et al., 2019). In contrast to Bereiter et al. (2018b) we do not impose a temporal change in this undersaturation, as we have no observational evidence about such a change in saturation state and we cannot provide a convincing argument on whether the saturation may have increased or decreased during glacial times. The overall slower ocean overturn in the glacial may suggest an increase in saturation, while the expansion of sea ice (especially in the Southern Ocean) would speak for a stronger undersaturation of the heavy noble gases. Accordingly, we refrain from changing the saturation state in our model. Sensitivity analyses show that a reduction of the glacial undersaturation to half its current value will lead to a warming of our glacial MOT reconstruction by a few tens of a degree.

Using this box model we can translate the reconstructed and corrected atmospheric elemental ratios ($\mathrm{Kr}/{\mathrm{N}}_{\mathrm{2}}$, $\mathrm{Xe}/{\mathrm{N}}_{\mathrm{2}}$, and $\mathrm{Xe}/\mathrm{Kr}$) into changes in MOT relative to the Holocene. Using our model-based correction approach described above leads to an MOT change over the last glacial–interglacial transition of on average 3.0 ± 0.4 ^∘C, 3.3 ± 0.4 ^∘C, and 3.6 ± 0.4 ^∘C for δKr $/$ N_2,atm, δXe $/$ N_2,atm, and $\mathit{\delta }\mathrm{Xe}/\mathrm{Kr}$_atm, respectively (see Fig. 6 and Table 2). This is in general agreement within uncertainties with independently derived results over the last glacial termination for the WAIS Divide ice core from Bereiter et al. (2018a).

Table 2 Absolute MOT anomaly based on $\mathrm{Kr}/{\mathrm{N}}_{\mathrm{2}}$ in degrees Celsius (^∘C) relative to the Holocene mean using both of our correction approaches. The Holocene value is the average of 10 EDC samples from 0–10.0 kyr BP, the LGM value is the average of nine samples between 17.8 and 26 kyr BP, LIG is the average of four samples between 120 and 129.9 kyr BP, and MIS 6 includes four samples from 135–160 kyr BP. MIS 9 and MIS 10 each comprise four samples.

As described above, the uncertainties in the model-based correction approach for elemental ratios are introduced by uncertainties in the modeled ΔT, the measurement errors of the elemental ratios, and the choice of ${\mathit{ϵ}}_{\text{Kr}/{\text{N}}_{\mathrm{2}}-{\text{N}}_{\mathrm{2}}}$. For the Ar isotopes we also corrected for the steady increase in ⁴⁰Ar through the decay of ⁴⁰K in the crust and mantle, which leads to a small long-term increase in ${\mathit{\delta }}^{\mathrm{40}/\mathrm{36}}\mathrm{Ar}$. This increase was experimentally determined by Bender et al. (2008) and has a slope of 0.066 ± 0.007 ‰ per million years, which we include in our error propagation. Moreover, the box model approach may add some minor systematic uncertainty, which we cannot quantify. Using the analytical uncertainties in the elemental ratios in the model-based approach leads to a stochastic uncertainty (1σ) in the MOT values between 0.35 and 0.41 ^∘C for all three elemental ratios and a systematic uncertainty of 0.38–0.51 ^∘C.

Figure 6 Reconstructed MOT for the EDC samples over the last 700 kyr relative to the Holocene value (using the model-based correction). MOTs are derived from $\mathit{\delta }\mathrm{Kr}/{\mathrm{N}}_{\mathrm{2}}$ (orange), $\mathit{\delta }\mathrm{Xe}/{\mathrm{N}}_{\mathrm{2}}$ (purple), and $\mathit{\delta }\mathrm{Xe}/\mathrm{Kr}$ (blue). The horizontal line in the middle of the filled box provides the corrected value. The boundaries of the filled boxes refer to the systematic uncertainty, while whiskers reflect the stochastic uncertainty introduced by the measurement errors (each 1σ). The red shaded intervals and the MIS numbers on top indicate interglacials as identified by Masson-Delmotte et al. (2010).

3 Post-coring gas loss effects on MOT reconstructions

Additional changes in the elemental ratios may occur if post-coring gas loss leads to fractionation between different gases. While visible cracks are avoided in the sample preparation, previous ice core studies showed that invisible microcracks can occur, leading to gas loss fractionation processes that affect, for example, $\mathit{\delta }{\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ and $\mathit{\delta }\mathrm{Ar}/{\mathrm{N}}_{\mathrm{2}}$, but little is known about gas loss effects on heavy noble gas ratios. Ice subject to such internal microcracks is most susceptible to post-coring gas loss, and accordingly we find systematic deviations in the MOT values derived from the three different elemental ratios in (i) ice from the bubble–clathrate transition zone (BCTZ), which is characterized by brittle ice, and (ii) in the oldest, bottommost meteoric ice of the EDC core, which is characterized by particularly large ice crystals and high in situ temperatures and which is prone to increased microfracturing. The following sections address potential gas loss fractionations in these ice zones.

3.1 Post-coring gas loss effects in the BCTZ

Post-coring fractionation of elemental gas ratios is caused by the selective loss of specific gases in poorly preserved samples or in samples from the bubble–clathrate transition zone (BCTZ), i.e., the depth section between the formation of the first clathrate and the disappearance of the last bubble (Lipenkov, 2000).

Figure 7 (a) Absolute MOT pair difference between two adjacent samples for the three different atmospheric noble gas ratios. (b) Difference in the MOT from $\mathit{\delta }\mathrm{Kr}/\mathrm{Xe}$ and $\mathit{\delta }\mathrm{Xe}/{\mathrm{N}}_{\mathrm{2}}$ compared to $\mathit{\delta }\mathrm{Kr}/{\mathrm{N}}_{\mathrm{2}}$. EDC, GRIP, NEEM, and Talos samples were measured at the University of Bern. WAIS samples were measured by Bereiter et al. (2018b). The grey shaded bars and arrows indicate the BCTZ for the different cores as identified by Neff (2014). All data are corrected for gravitational and thermal fractionation using argon isotopes. EDC samples are also corrected for HID fractionation. (c) Potential xenon change required for the EDC samples to match the MOTs using a least square approach. Given in red is the potential xenon change that minimizes the difference between the MOTs derived from $\mathrm{Xe}/\mathrm{Kr}$, $\mathrm{Kr}/{\mathrm{N}}_{\mathrm{2}}$, and $\mathrm{Xe}/{\mathrm{N}}_{\mathrm{2}}$. Negative values indicate a xenon loss, i.e., a depletion of the raw elemental ratios $\mathrm{Xe}/{\mathrm{N}}_{\mathrm{2}}$ and $\mathrm{Xe}/\mathrm{Kr}$ caused by post-coring gas loss. Positive values represent apparent net xenon “gains” due to higher losses of krypton and nitrogen compared to xenon. (d) Evolution of the grain radius in the EDC core (EPICA community members, 2004; Durand et al., 2009). Note that the ages provided on the upper x axis only apply for samples from the EDC ice core. MIS interglacials (red shaded intervals) as identified by Masson-Delmotte et al. (2010) are labeled at the top for orientation.

Bender (2002) called such gas loss through cracks or microcracks “core-cracking fractionation” and observed this effect to strongly affect ${\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ in the BCTZ of the Vostok ice core. The effect is explained by gas fractionation occurring between coexisting bubbles and clathrates in combination with post-coring gas loss. In our record we observe small but linearly increasing differences with age (and depth) in the MOT derived from the different noble gas combinations in ice at the top of the BCTZ, i.e., in the zone where clathrates form (600–775 m; 24–40 kyr BP). In particular, MOTs derived from $\mathrm{Kr}/{\mathrm{N}}_{\mathrm{2}}$ show higher temperatures compared to $\mathrm{Xe}/{\mathrm{N}}_{\mathrm{2}}$ and $\mathrm{Xe}/\mathrm{Kr}$ in this ice, while the three values agree within their uncertainties in pure bubble ice. Similar systematic offsets are observed around the BCTZ in cores from various sites (GRIP, NEEM, and Talos) but to a variable extent, as shown in Fig. 7b. Bereiter et al. (2018a) also observed pronounced alteration in the trapped gas towards the lower end of the BCTZ or even below for the WAIS ice core.

In the following discussion, we will try to motivate the processes that we think are responsible for this fractionation. Studies of air bubbles and clathrate inclusions on the Dome Fuji ice core provide information on clathrate formation and clathrate growth in cold Antarctic ice (Uchida et al., 2011). Clathrate formation starts when the bubble pressure, which increases linearly with depth due to the increasing hydrostatic pressure, reaches the dissociation pressure of air hydrate. The dissociation pressure of the air hydrate is given by the dissociation pressure of hydrates formed by the pure gas species multiplied by the mole fraction of the gas species in air (Miller 1969). Accordingly, the air hydrate dissociation pressure in the ice is mainly determined by the dissociation pressure of N₂ and O₂, with minor contributions from Ar or CO₂. We assume that any SII-clathrate-forming gas species (Sloan and Koh, 2007) will be quantitatively incorporated into the SII air hydrate structure during its nucleation from an air bubble; however, we speculate that SI-forming gases (such as Xe) may not (or not quantitatively) be incorporated into the SII air hydrate structure in the ice.

As described in Uchida et al. (2011), the number and also the radius of air hydrate crystals increase in the BCTZ at Dome Fuji, which has very similar climatic conditions compared to EDC. This can be interpreted in the BCTZ as more and more hydrates slowly nucleating from bubbles and bubbles and clathrates coexisting for a long time. Due to the higher partial pressures of gas species in bubbles compared to clathrates a constant permeation flux of gas from the bubbles to the clathrates exists, which is the main cause of clathrate growth in the BCTZ at Dome Fuji (or EDC). Note that due to the different permeation constants of different gas species in ice, this transport implies a fractionation of different gas species, which enriches fast-permeating gases (such as O₂; Ikeda et al., 2000) in the hydrate and depletes them in the remaining bubble. For Xe we speculate that there is either no permeation flux (if Xe is not included in the SII hydrate structure) or the permeation flux is very small due to the very low permeation coefficient of the large Xe atom. Thus, Xe is depleted in the air hydrate forming in the BCTZ and enriched in the remaining bubbles. After complete collapse of the bubble, we expect the Xe to be “dissolved” in the ice matrix.

Below the BCTZ Uchida et al. (2011) observe an increase in hydrate radius and at the same time a decline in hydrate number. This can be explained by an Ostwald ripening process, leading to the growth of large hydrates at the cost of small ones through the Gibbs–Thomson effect (Uchida et al., 2011). Again, the accompanying gas transport from a small to a large hydrate crystal is supplied by (fractionating) permeation of gas species through the ice. If this process goes on for a sufficiently long time, however, a diffusive re-equilibration of the mole fractions of gas species in the hydrates (such as CO₂) is been observed (Lüthi et al., 2010). We assume that this permeation acts on all SII-forming gas species but that SI-forming Xe either does not permeate at all or that the permeation flux is strongly reduced compared to other gas species. Thus, Xe would be depleted in the large clathrates formed by this Ostwald ripening process and would remain dissolved in the ice matrix after complete disappearance of a small clathrate.

We attribute the gradually increasing offsets in reconstructed MOT at the top of the BCTZ to core-cracking gas loss in combination with the fractionation of gases by permeation. Gas enclosed in bubbles is more likely to be lost through cracks after core retrieval than gas in hydrates; thus, if post-coring gas loss from bubbles occurs, the remaining gas of the ice sample will be enriched in O₂ relative to N₂ and depleted in Xe. The small but significant $\mathit{\delta }\mathrm{Xe}/\mathrm{Kr}$ and $\mathit{\delta }\mathrm{Xe}/{\mathrm{N}}_{\mathrm{2}}$ decrease in the BCTZ can therefore be explained by the observation that during the gas loss process, Xe is preferentially lost via cracks from the Xe-enriched bubbles.

To quantify a potential Xe loss we use a least squares approach that artificially changes the Xe abundance in the ice to minimize the difference between the three MOT proxies. The respective calculated post-coring Xe loss is shown in Fig. 7c, where the per mill value is the change in the elemental ratio used to minimize the MOT differences compared to $\mathrm{Kr}/{\mathrm{N}}_{\mathrm{2}}$. It shows that Xe loss is almost negligible in the top 500 m, i.e., in bubble ice, and increases to about 2 ‰, meaning that the $\mathrm{Xe}/{\mathrm{N}}_{\mathrm{2}}$ elemental ratio is depleted in Xe by up to 2 ‰ caused by post-coring Xe loss from ice in the BCTZ.

Previous studies have shown that the gas loss effect on ${\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ is most pronounced when the ice is brittle with an enhanced number of cracks but that the quality of the data is also affected in non-brittle ice if the samples are stored for an extended amount of time at relatively warm temperatures of around −25 ^∘C (Ikeda-Fukazawa et al., 2005; Landais et al., 2012). The storage temperature may also affect our noble gas ratio results because our samples from the last transition were stored for about 15 years at about −22 ^∘C. The samples from older climatic intervals presented in this study were stored at Dome C and later in Bern in a freezer at a temperature below −50 ^∘C. Thus, no gas loss is expected for these samples during storage. However, these samples experienced a short-term warming during transport, when cooling failed during transit. Temperatures were as warm as −16 ^∘C at the bottom of the well-insulated transport box and −6 ^∘C on top of the box for several hours; however, the ice never reached the melting point.

Figure 8 Measurements of δ${\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ and $\mathit{\delta }\mathrm{Ar}/{\mathrm{N}}_{\mathrm{2}}$ from the Dome C ice core on the EDC depth scale. The δ${\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ record from Extier et al. (2018) is corrected for gravitational enrichment using δ¹⁵N₂. Our δ${\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ data are corrected for gravitational enrichment using ${\mathit{\delta }}^{\mathrm{40}/\mathrm{36}}\mathrm{Ar}$, while the $\mathit{\delta }\mathrm{Ar}/{\mathrm{N}}_{\mathrm{2}}$ data are corrected for gravitational enrichment, thermal diffusion, kinetic fractionation, and geological Ar outgassing. The data show a large variability in the BCTZ and reproducible values for clathrated ice. (a) The green line represents measurements from Extier et al. (2018), and diamonds represent the EDC measurements of δ${\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ in this study. (b) EDC measurements of $\mathit{\delta }\mathrm{Ar}/{\mathrm{N}}_{\mathrm{2}}$ from this study. Interglacial MISs (red shaded areas) as identified by Masson-Delmotte et al. (2010) are labeled for orientation.

To assess the potential loss of gas due to high storage or transport temperature, we compare the ${\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ that we measure alongside the noble gases with previous ${\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ measurements on EDC ice stored at −50 ^∘C by Extier et al. (2018); the latter did not experience significant warming during transport. Similar to previous observations in BCTZ ice, our ${\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ values suggest a variable enrichment of up to 4 ‰ in O₂ relative to N₂ at the top of the BCTZ and also slightly elevated $\mathrm{Ar}/{\mathrm{N}}_{\mathrm{2}}$ relative to the values in pure bubble ice (as shown in Fig. 8). This suggests that our samples from the BCTZ are indeed subject to slight core-cracking fractionation. Note that for $\mathrm{Kr}/{\mathrm{N}}_{\mathrm{2}}$ we do not see such an enrichment at the top of the BCTZ, which implies that the difference in permeation constants for Kr and N₂ is too small to create a sizable anomaly. For the deeper ice no significant offset relative to the ${\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ data from Extier et al. (2018) can be found, suggesting that the short cooling interruption during transport did not affect our elemental ratios. For the $\mathrm{Ar}/{\mathrm{N}}_{\mathrm{2}}$ ratios of samples older than 40 kyr BP we see a number of glacial values that are more depleted than the average value during the Last Glacial Maximum (LGM); however, the large spread of the data does not allow us to draw a firm conclusion on an Ar loss for these older samples. Furthermore, both ${\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ and $\mathrm{Ar}/{\mathrm{N}}_{\mathrm{2}}$ of neighboring samples agree well with each other, indicating a common cause for the large overall variability. The main factor influencing ${\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ on these timescales is a change in local summer insolation (Kawamura et al., 2007; Extier et al., 2018), and $\mathrm{Ar}/{\mathrm{N}}_{\mathrm{2}}$ should be affected in a similar fashion. This suggests that gas loss during transport and storage is of minor importance, as such artifactual gas loss is typically a chaotic process that randomly affects some samples more than others. Based on this, we refrain from making a gas loss correction for ${\mathit{\delta }}^{\mathrm{40}/\mathrm{36}}\mathrm{Ar}$ as described in Sect. 2.2.6.

In summary, although the Xe depletion is relatively small in the ice at the top of the BCTZ, we believe that the MOTs derived from our $\mathrm{Xe}/\mathrm{Kr}$ and $\mathrm{Xe}/{\mathrm{N}}_{\mathrm{2}}$ values are slightly biased towards overly cold temperatures. Accordingly, we use the MOT values derived from $\mathrm{Kr}/{\mathrm{N}}_{\mathrm{2}}$ to quantitatively reconstruct changes in MOT for ice from the BCTZ, as $\mathrm{Kr}/{\mathrm{N}}_{\mathrm{2}}$ appears not to be affected by gas loss and as $\mathrm{Kr}/{\mathrm{N}}_{\mathrm{2}}$ replicates agree well within their analytical uncertainty thus making an effect of (unreproducible) gas loss unlikely. The latter also suggests that any fractionation of Kr relative to N₂ that may be caused by a permeation fractionation between bubbles and clathrates is too small to have a sizable effect on our measured values.

3.2 Post-coring gas loss in ice deeper than 2500 m

In clathrate ice below the BCTZ, MOT values derived from the three different elemental ratios showed good agreement down to an age of about 450 kyr BP but revealed systematic offsets mainly for the lowest 400 m, where the MOTs derived from $\mathit{\delta }\mathrm{Xe}/\mathrm{Kr}$ and $\mathit{\delta }\mathrm{Xe}/{\mathrm{N}}_{\mathrm{2}}$ suggest much higher temperatures than those derived from $\mathit{\delta }\mathrm{Kr}/{\mathrm{N}}_{\mathrm{2}}$. A similar observation can be made for two interglacial samples from MIS 9. Moreover, the values derived from $\mathit{\delta }\mathrm{Xe}/\mathrm{Kr}$ and $\mathit{\delta }\mathrm{Xe}/{\mathrm{N}}_{\mathrm{2}}$ of neighboring quasi-replicate samples in the lowest 400 m of the core do not agree within their analytical uncertainties, while those for $\mathit{\delta }\mathrm{Kr}/{\mathrm{N}}_{\mathrm{2}}$ do, as shown by differences in the MOT of each sample pair in Fig. 7a. Each of the $\mathit{\delta }\mathrm{Kr}/{\mathrm{N}}_{\mathrm{2}}$ pairs (with the exception of one sample at a depth of about 3100 m from 630 kyr BP) agrees within less than 0.5 ^∘C, in line with our stochastic error estimate (see Fig. 7). The cause of the outlier in $\mathit{\delta }\mathrm{Kr}/{\mathrm{N}}_{\mathrm{2}}$ at 630 kyr BP has not yet been identified. Interestingly, $\mathit{\delta }\mathrm{Kr}/{\mathrm{N}}_{\mathrm{2}}$ pair differences in the bottommost ice are also not elevated compared to shallower fully clathrated ice samples in the time interval 100–450 kyr BP. Taking the good reproducibility of $\mathit{\delta }\mathrm{Kr}/{\mathrm{N}}_{\mathrm{2}}$ and the elevated MOTs derived from $\mathit{\delta }\mathrm{Xe}/\mathrm{Kr}$ and $\mathit{\delta }\mathrm{Xe}/{\mathrm{N}}_{\mathrm{2}}$ in the deepest ice together, we conclude that a Kr and N₂ core-cracking loss must have occurred, which leaves the ice enriched in Xe.

In the EDC core, Durand et al. (2009) observed large variations in the fabric of the ice in different samples and extraordinarily large grains in interglacial ice below 2500 m, as shown in Fig. 7. A recent microscopic inspection of two deep EDC ice samples stored at the Alfred Wegener Institute in Bremerhaven at −30 ^∘C showed significant visual glittering in this ice, indicative of so-called plate-like inclusions (PLIs) due to clathrate relaxation (Bereiter et al., 2015). These suggest the coexistence of reformed gaseous inclusions and clathrates in this increasingly relaxed ice. Weikusat et al. (2012) and Nedelcu et al. (2009) studied the ice structure and relaxation of the EDML core. Although the ice is comparatively young in this core and the accumulation rates are larger than at the EDC site, it is likely that the observations of EDML also apply to some degree to the EDC core. Particularly interesting are relaxation-induced air inclusions such as microscopic bubbles (“microbubbles”) and PLIs that are found with increasing depth at grain boundaries and clathrate hydrate surfaces (Nedelcu et al., 2009). The PLIs are flat inclusions of air (a few micrometer in thickness) and are more common in deep, fully clathrated ice, while the microbubbles are spherical objects and more common in shallower, bubbly ice (Weikusat et al., 2012). Direct measurements of such microscopic air inclusions are only possible for the main air constituents N₂ and O₂ using Raman spectroscopy. The results from Nedelcu et al. (2009) show substantially higher ${\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ values in such inclusions compared to air, indicating an enrichment of O₂ in these gas-filled features similar in magnitude to the O₂ enrichment reported for air hydrates in the BCTZ of polar ice cores.

Diffusive gas loss from the outer surface of our samples can be safely ruled out due to the large sample size and the careful removal of several millimeters of ice from each surface (Bereiter et al., 2009). In combination with the complete gas extraction achieved in our melt extraction technique, any occurring net gas fractionation in our samples requires a gas loss process through microcracks and some gas fractionating exchange between gas enclosures in the ice. Similar to the core-cracking gas loss in the BCTZ discussed above, any gas loss through microcracks in the ice is more likely for gas-filled features such as bubbles and PLIs than for clathrates. Hence, we would assume that gas loss from O₂-enriched PLIs would lead to a depletion in the ${\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ values measured in our ice samples. As shown in Fig. 8, however, an O₂ enrichment is observed that is increasing for older ice and that has been attributed to a long-term change in the O₂ content of the atmosphere over the last 800 kyr (Stolper et al., 2016). Thus, either an O₂ loss did not occur or it is very small and dwarfed by the atmospheric O₂ change.

Nevertheless, we suggest that the apparent enrichment of Xe relative to N₂ and Kr in our data is related to a fractionation and gas loss process that is opposite to that what is observed in the BCTZ and that is related to the clathrate relaxation process. The fractionation is apparently particularly pronounced in ice below 2500 m with large crystals (see Fig. 7d) and in situ temperatures warmer than −10 ^∘C, at which this ice may experience enhanced microcracking during core retrieval and storage. Once the ice core is brought to the surface, it starts to relax and gas-filled features such as PLIs slowly form (the higher the storage temperature of the core, the faster the relaxation). Note that due to the very cold storage temperature of our samples older than 40 kyr ($<-$50 ^∘C) at Concordia Station and later in a −50 ^∘C freezer chest in our lab, relaxation of the core is strongly suppressed, and we expect far fewer PLIs to have formed in our samples compared to the EDML samples in Nedelcu et al. (2009). Thus, it is likely that any O₂ loss is too small to be discerned.

How can we explain a Xe enrichment in our samples without significant ${\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$ fractionation? As discussed in Sect. 3.1, the SII clathrates are likely to be depleted in Xe and part of the Xe is expected to remain dissolved in the ice matrix. Thus, while any gas inclusions reforming from clathrates will be depleted in Xe and can be subject to gas loss, the Xe molecules dissolved in the ice are not accessible to gas loss. Accordingly, we speculate that any gas loss leads to an enrichment of Xe in our deep EDC ice samples, while all other gas species do not show significant fractionation as long as very cold storage ensures the clathrate relaxation to be small. This is also the reason why $\mathrm{Kr}/{\mathrm{N}}_{\mathrm{2}}$ ratios apparently do not show any signs of gas loss fractionation as illustrated by the small pair differences between adjacent samples in the deep ice. Experiments using old EDC ice that has been stored for an extended time at a higher storage temperature may help researchers to test this hypothesis in the future.

In summary, apart from analytical issues related to achieving the required high precision and accuracy in atmospheric values of elemental ratios, fractionation and gas loss effects also have to be taken into account when interpreting MOT values in deep ice cores. As observed earlier for δ${\mathrm{O}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}$, the gas composition in brittle ice from the BCTZ suffers from post-coring gas loss effects, while pure bubble ice provides good storage conditions for permanent gases. MOT values from pure clathrate ice below the BCTZ but well above the bedrock (>400 m) show little signs of post-coring gas loss and provide reliable MOTs for all three elemental ratios. In the warm deep ice closer to the bedrock, we observed an apparent Xe enrichment, and the $\mathrm{Xe}/\mathrm{Kr}$ and $\mathrm{Xe}/{\mathrm{N}}_{\mathrm{2}}$ ratios did not allow us to derive reliable MOT values. The consistency of $\mathit{\delta }\mathrm{Kr}/{\mathrm{N}}_{\mathrm{2}}$ values over the entire record, particularly their reproducibility within the analytical uncertainty, gives us confidence that reliable MOT values can be derived from this parameter throughout the EDC ice core.

4 MOT snapshots of peak glacial and interglacial periods during the last 700 kyr

Following the discussion about gas loss above, MOTs derived from $\mathit{\delta }\mathrm{Xe}/{\mathrm{N}}_{\mathrm{2}}$ and $\mathit{\delta }\mathrm{Xe}/\mathrm{Kr}$ show significant anomalies in samples from the BCTZ, especially from very deep ice as shown in Fig. 6.

Based on the good reproducibility of $\mathit{\delta }\mathrm{Kr}/{\mathrm{N}}_{\mathrm{2}}$ values within the analytical uncertainty throughout the record, we conclude that the MOT derived from $\mathit{\delta }\mathrm{Kr}/{\mathrm{N}}_{\mathrm{2}}$ is not significantly affected by the gas loss described in Sect. 3 and shows consistent results throughout the record. One sample at 638 kyr BP is, however, considered an outlier for two reasons: the sample has a pair difference in $\mathit{\delta }\mathrm{Kr}/{\mathrm{N}}_{\mathrm{2}}$, which is larger than 3σ of the reproducibility, and the order of the three proxies indicates a substantial xenon loss in contrast to all other samples from pure clathrate ice. We neglect the value at 638 kyr BP and restrict our discussion of past MOT changes to the values derived from $\mathit{\delta }\mathrm{Kr}/{\mathrm{N}}_{\mathrm{2}}$ as presented in Fig. 9.

Figure 9 MOT record relative to the Holocene from the $\mathrm{Kr}/{\mathrm{N}}_{\mathrm{2}}$ ratios in comparison with other climate records. The EDC record is based on 85 samples, which are displayed as orange boxes (except one outlier at 638 kyr BP indicated in grey). The uncertainty is shown as described in Sect. 2.3. (a) Deep-ocean reconstructions from a global compilation of 49 paired planktic δ¹⁸O records are shown in dark blue (Shakun et al., 2015). Given in cyan is the benthic $\mathrm{Mg}/\mathrm{Ca}$-derived temperature from Elderfield et al. (2012). (b) The CO₂ record is based on measurements from the EDC and Vostok ice core (Nehrbass-Ahles et al., 2020; Bereiter et al., 2015; Lüthi et al., 2008; Petit et al., 1999). (c) Antarctic temperatures are taken from Parrenin et al. (2013). Red shaded bars and MIS numbers on top indicate interglacials according to Masson-Delmotte et al. (2010).

Despite the pointwise character of our data, the $\mathit{\delta }\mathrm{Kr}/{\mathrm{N}}_{\mathrm{2}}$-derived MOT record generally supports the main findings from the deep-ocean temperature reconstruction of Elderfield et al. (2012) and Shakun et al. (2015), but in contrast to these studies our record is based on a purely physical proxy and by definition is globally representative for MOT. The main findings are the following.

• i.

  Some episodes during the interglacial periods MIS 5.5 and MIS 9.3 are significantly warmer than the Holocene.

• ii.

  Glacial MOTs are similar throughout the last 800 kyr, although glacials prior to 450 kyr BP appear slightly colder.

• iii.

  The MOTs during interglacial periods prior to 450 kyr BP are significantly colder than after the Mid-Brunhes Event.

Note that each of our samples only represents an MOT snapshot representative of a time interval of a few centuries (Fourteau et al., 2019) given by the width of the gas age distribution for the EDC ice core. Thus individual samples are not representative of the mean interglacial (or glacial) value but may be subject to the millennial-scale variability of MOT. We stress that MOT is strongly influenced by changes in the Atlantic Meridional Overturning Circulation (AMOC) (Baggenstos et al., 2019; Shackleton et al., 2020) and can therefore vary on multi-centennial to millennial timescales, as also seen in coupled climate model experiments (Pedro et al., 2018; Galbraith et al., 2016). A sudden slowdown of the AMOC leads to a long-lasting accumulation of heat in the ocean interior (Pedro et al., 2018; Galbraith et al., 2016; Ritz et al., 2011), while a recovery of the AMOC slowly lowers the MOT. Accordingly, samples very close to the onset of interglacials, which may still be affected by a previously reduced AMOC, may still show elevated MOTs. The heat stored in the global ocean is only removed in the course of 2000–3000 years into the subsequent interglacial after strengthening of the AMOC. This is clearly observed in higher-resolution data for MIS5.5 (Shackleton et al., 2020), which allowed us to temporally resolve this transient feature. Such a deepwater mass reorganization that marks the onset of the subsequent interglacial period is also inferred from changing ocean water tracer distributions (Deaney et al., 2017).

In line with this, our MOT data based on $\mathit{\delta }\mathrm{Kr}/{\mathrm{N}}_{\mathrm{2}}$ from the onset of MIS 5.5 are on average 1.0 ± 0.4 ^∘C warmer than the Holocene and decline to interglacial values only 0.7 ^∘C warmer than the Holocene later in MIS5.5. This is in agreement with the findings from Shackleton et al. (2020) that show significantly warmer MOTs reaching a maximum value of 1.1 ± 0.3 ^∘C at about 129 kyr BP. Shackleton et al. (2020) attribute the early last interglacial maximum in MOT to the weakening of the AMOC that occurred over termination II during the major iceberg discharge event from Hudson Strait in the North Atlantic commonly known as Heinrich event 11 (Capron et al., 2014). The subsequent release of heat from the ocean after the AMOC resumed led to a MOT decrease of ∼ 1 ^∘C with stable MOT only reached after 127 kyr BP.

A similar but even more pronounced ocean warming feature is found in our new data during interglacial MIS 9.3. The samples later in the interglacial (329 and 330 kyr BP) show within uncertainty the same temperatures as the Holocene; however, the very early interglacial samples (334 and 335 kyr BP) show 1.8 ± 0.4 ^∘C warmer MOTs than the Holocene. Similarly, global deep-ocean temperatures for early MIS 9.3 compiled from marine sediment records in the study by Shakun et al. (2015) are intermittently nearly 2 ^∘C warmer than the Holocene. Our two early MIS 9.3 samples occur also in parallel with a rapid CO₂ and CH₄ increase at the very end of termination IV. At the onset of interglacial MIS 9.3, Nehrbass-Ahles et al. (2020) measured the highest CO₂ values during the last 800 kyr, coeval with the highest CH₄ levels (Loulergue et al., 2008) and warm Antarctic temperatures, providing an excellent chronological constraint for the timing of AMOC resumption. An intriguing hypothesis for the high values at the onset of interglacial MIS 9.3 is therefore that the MOT overshoot observed in the final phases of termination IV is also a result of a long preceding period of freshwater-induced AMOC suppression and deepwater mass reorganization (Ganopolski and Brovkin, 2017) that only ended with the AMOC resumption at that time. This suggests, similar to the case in MIS 5.5, that MOT was still elevated by the reduced AMOC in the millennia before and was able to regain interglacial levels several millennia after the rapid CO₂ and CH₄ increases. This emphasizes the transient character of this feature and the role of an increasing AMOC in Antarctic and ocean temperature as well as in marine and terrestrial biogeochemical cycles.

Our data show glacial mean ocean temperatures very similar throughout the last 800 kyr. On average our measured glacial MOTs are 3.3 ± 0.4 ^∘C colder than the Holocene. This is in agreement with the deep-ocean anomaly derived from Elderfield et al. (2012) based on a sediment record from Chatham rise east of New Zealand. There appears to be a slight increasing trend in glacial MOT over the last 700 kyr that is also suggested by the global ocean data compilation of Shakun et al. (2015); however, the trend in our data is of the same order as the uncertainty of our reconstruction. Accordingly, we refrain from interpreting this trend at this point in time; more higher-resolution MOT data are required to show whether this feature persists.

Our data also support the finding of lukewarm interglacials in MOT prior to the Mid-Brunhes Event 450 kyr ago (Shakun et al., 2015; Elderfield et al., 2012), as also seen in Antarctic temperatures (EPICA community members, 2004), but also lower interglacial concentrations of the greenhouse gases CO₂ (Bereiter et al., 2015; Lüthi et al., 2008) and CH₄ (Loulergue et al., 2008), at least for MIS 13–17. Thus, a potential compensation for the excess effective energy not used for ice sheet melting during lukewarm interglacials, as discussed in the Introduction, can be ruled out based on our ice core and the marine sediment data. Accordingly, the lukewarm interglacials must be connected to a net change in the radiative balance of the globe compared to later warmer interglacials. The MOT anomaly relative to the Holocene during the interglacials MIS 13, 15.1, 15.5, and 17 is −1.6 ± 0.4, −1.7 ± 0.4, −1.3 ± 0.4, and −1.8 ± 0.4 ^∘C, respectively. Interglacial bottom-water temperatures also increased by about 1 ^∘C after ∼ 450 kyr BP (Elderfield et al., 2012; Shakun et al., 2015); however, Holocene bottom-water temperatures seem to be colder again than during the three preceding interglacials MIS 5.5, MIS 9.3, and MIS 11. Using the global benthic δ¹⁸O stack by Lisiecki and Raymo (2005), Bintanja et al. (2005) reconstructed global ice volume changes and the contributions of individual ice sheets back in time. Again, the lukewarm interglacials are significantly different compared to the later interglacials in this model reconstruction and are characterized by up to 20 m lower sea level compared to the Holocene, with the largest contribution to this sea level reduction caused by larger ice sheets in Eurasia. Note, however, that sea levels lower than in the Holocene during lukewarm interglacials can be observed in some but not all sea level reconstructions from marine sediments (Rohling et al., 2014; Spratt and Lisiecki, 2016). In any case, an increase in ice volume during lukewarm interglacials cannot be easily attributed to changes in orbital forcing, as the effective energy at high northern latitudes (Tzedakis et al., 2017) during these intervals did not significantly differ from the insolation conditions during later interglacials. Thus, the lukewarm interglacials require additional changes in the Earth system that prevented Northern Hemisphere ice sheets from retreating and temperatures from increasing to the same extent as in later interglacials.

Figure 10 Scatter plot showing a linear relationship between logarithmic CO₂ (values taken from Bereiter et al., 2015) and the MOT relative to the Holocene. MOT reconstructions are based on $\mathrm{Kr}/{\mathrm{N}}_{\mathrm{2}}$ and corrected for thermal, gravitational, and HID fractionation. MOTs from the last 25 kyr are shown in blue, samples from 25–40 kyr BP are shown in orange, and samples from 120 to 700 kyr BP are in green. The linear regression (black line) is based on EDC data only and has an r² of 0.75. Given in grey is the MOT reconstruction over the last 22 kyr from Bereiter et al. (2018b) relative to the Holocene.

Irrespective of the ultimate cause of the incomplete deglaciation during the lukewarm interglacials, our significantly lower MOT values during these lukewarm interglacials provide robust evidence for a significantly altered energy budget of our planet compared to the Holocene. Note that an MOT lowering of about 1 ^∘C translates into a decrease in the ocean heat content of 5–6×10²⁴ J, or about 30 % of our typical glacial–interglacial change in heat content. We can estimate whether the deviation in climate boundary conditions (greenhouse gas concentrations and albedo) is sufficient to explain the heat deficit. Lukewarm interglacial CO₂ concentrations are typically 245 ppm (instead of 280 ppm in later interglacials) and CH₄ concentrations 620 ppb (instead of 700 ppb thereafter). This translates into a radiative forcing of about −0.8 W m⁻² (Myhre et al., 1998), and using a recent paleo-based value of the equilibrium climate sensitivity of 3.8 ^∘C for a doubling of CO₂ (Tierney et al., 2020) (which includes fast feedbacks such as sea ice coverage) it translates into a reduction of the global mean surface temperature (GMST) of about 0.8 ^∘C. In addition, we also need to take into account the long-term Earth system feedbacks such as planetary albedo increase due to extended land ice coverage. For the LGM, the fraction of greenhouse gas forcing to total forcing has been estimated as approximately 0.4 (Baggenstos et al., 2019). Using the greenhouse gas forcing for the lukewarm interglacials as given above, the estimated total forcing is thus close to 2 W m⁻², which yields an expected GMST decrease of 2.0 ^∘C for the lukewarm interglacials. However, comparing global surface temperature to MOT is not straightforward (Bereiter et al., 2018b) because (i) MOT is biased towards sea surface temperatures in deepwater formation regions, for example in the North Atlantic, and is thus subject to polar amplification, and (ii) land temperature change is poorly known and cannot easily be linked to MOT change. However, a recent study has tightly constrained LGM GMST change to 6 ^∘C (Tierney et al., 2020), which goes along with an approximately 3 ^∘C MOT change for the LGM. Assuming this relationship also holds for earlier interglacials, we scale our average MOT decrease of 1.6 ^∘C into a GMST reduction of 3.2 ^∘C. This is somewhat larger than the expected GMST lowering of 2.0 ^∘C due to radiative forcing and suggests that either the ratio of greenhouse gas forcing to total forcing or the scaling from MOT to GMST was different during the lukewarm interglacials. Note, however, that both of these factors are poorly constrained, and, for example, using 0.3 instead of 0.4 for the forcing ratio along with a slightly smaller MOT and GMST scaling would bring the two estimates into agreement. Thus, our MOT lowering during lukewarm interglacials is within the range expected from the radiative forcing caused by the lower CO₂ and CH₄ in those intervals and the connected longer-term Earth system feedbacks.

The close relationship of MOT and global radiative forcing with changing greenhouse gas concentrations is not only valid for lukewarm vs. full interglacials but to first order also holds for the entire noble-gas-based MOT record available so far. Bereiter et al. (2018b) showed that MOT is closely correlated with Antarctic temperature and the MOT evolution is essentially synchronous with CO₂ during the last transition. Our record from the last 40 kyr supports this correlation. The scatter plot presented in Fig. 10 shows the correlation between MOT and logarithmic CO₂, as expected from the radiative forcing of surface temperatures by CO₂ concentration (IPCC, 2013). The squared correlation coefficient for all the reconstructed MOTs from EDC and logarithmic CO₂ suggests that 75 % of the variance in the MOT record can be explained by changing CO₂ and the accompanying Earth system feedbacks. Taking a closer look at Fig. 10 we see again that the peak glacials and the lukewarm interglacials prior to 450 kyr BP show approximately 1 ^∘C colder temperatures than expected from the CO₂ relationship derived for the last 40 kyr. In view of the considerable systematic uncertainties introduced by the correction of the data, we refrain from speculating on a possible cause of these lower temperatures at this point.

5 Conclusions

Our study presents MOT snapshots of peak glacial and interglacial periods during the last 700 kyr. The comparison of MOT values derived from different elemental ratios revealed systematically lower MOTs based on $\mathit{\delta }\mathrm{Xe}/{\mathrm{N}}_{\mathrm{2}}$ and $\mathit{\delta }\mathrm{Xe}/\mathrm{Kr}$ values from the BCTZ and higher MOTs from $\mathit{\delta }\mathrm{Xe}/{\mathrm{N}}_{\mathrm{2}}$ and $\mathit{\delta }\mathrm{Xe}/\mathrm{Kr}$ in the ice deeper than 2500 m. We attribute this to a Xe depletion during clathrate formation through its low permeation coefficient and/or reduced incorporation of Xe in SII clathrates in the ice. This leads to an enrichment of Xe in bubbles in the BCTZ and a depletion of Xe in reforming PLIs in the deep ice. Together with gas loss through post-coring microcracks, this leads to a depletion of Xe relative to N₂ and Kr in the BCTZ and an enrichment in the deep ice. Continuing work on understanding firn gas, especially the HID, is crucial to reduce the uncertainties of the MOT reconstructions based on noble gases.

Despite these limitations we find that glacial MOTs based on $\mathit{\delta }\mathrm{Kr}/{\mathrm{N}}_{\mathrm{2}}$, which appear to be free of gas loss and fractionation processes, are on average 3.3 ± 0.4 ^∘C colder than the Holocene throughout this noble gas record and that the MOTs during interglacial periods prior to 450 kyr BP are on average 1.6 ± 0.4 ^∘C lower compared to the Holocene. Episodes during the onsets of interglacials MIS 5.5 and 9.3 are significantly warmer than the Holocene and can be attributed to the ocean heat uptake during preceding times of lowered AMOC. The most pronounced feature of the record is the significantly increased MOT at the onset of MIS 9.3, which is coeval with the CO₂ and CH₄ overshoot at that time and reflects maximum ocean heat uptake before the AMOC resumption.

Our data provide only snapshots of MOT in glacial and interglacial conditions; however, these define the natural range in which the global heat budget has changed over the last 700 kyr. Our data show that the range in MOT can be quantitatively explained by changes in global radiative forcing (due to changes in greenhouse gas concentrations), taking into account long-term Earth system feedbacks, and that MOT is biased towards sea surface temperature changes in high latitudes. More high-resolution and high-precision MOT data from noble gases in ice cores (avoiding gas loss issues) and a better understanding of the firn fractionation processes affecting the gas composition in ice cores are needed to make use of the full potential of noble-gas-derived MOT in future studies.