Articles | Volume 16, issue 7
The Cryosphere, 16, 2947–2966, 2022
The Cryosphere, 16, 2947–2966, 2022
Research article
22 Jul 2022
Research article | 22 Jul 2022

Gas isotope thermometry in the South Pole and Dome Fuji ice cores provides evidence for seasonal rectification of ice core gas records

Gas isotope thermometry in the South Pole and Dome Fuji ice cores provides evidence for seasonal rectification of ice core gas records
Jacob D. Morgan1, Christo Buizert2, Tyler J. Fudge3, Kenji Kawamura4,5,6, Jeffrey P. Severinghaus1, and Cathy M. Trudinger7 Jacob D. Morgan et al.
  • 1Scripps Institution of Oceanography, University of California, San Diego, La Jolla, CA 92093, USA
  • 2College of Earth, Ocean, and Atmospheric Sciences, Oregon State University, Corvallis, OR 97331, USA
  • 3Department of Earth and Space Sciences, University of Washington, Seattle, WA 98195, USA
  • 4National Institute of Polar Research, Tokyo 190-8518, Japan
  • 5Department of Polar Science, The Graduate University of Advanced Studies (SOKENDAI), Tokyo 190-8518, Japan
  • 6Japan Agency for Marine-Earth Science and Technology (JAMSTEC), Yokosuka 237-0061, Japan
  • 7Climate Science Centre, CSIRO Oceans and Atmosphere, Aspendale, Victoria 3195, Australia

Correspondence: Jacob D. Morgan (


Gas isotope thermometry using the isotopes of molecular nitrogen and argon has been used extensively to reconstruct past surface temperature change from Greenland ice cores. The gas isotope ratios δ15N and δ40Ar in the ice core are each set by the amount of gravitational and thermal fractionation in the firn. The gravitational component of fractionation is proportional to the firn thickness, and the thermal component is proportional to the temperature difference between the top and bottom of the firn column, which can be related to surface temperature change. Compared to Greenland, Antarctic climate change is typically more gradual and smaller in magnitude, which results in smaller thermal fractionation signals that are harder to detect. This has hampered application of gas isotope thermometry to Antarctic ice cores.

Here, we present an analytical method for measuring δ15N and δ40Ar with a precision of 0.002 ‰ per atomic mass unit, a two-fold improvement on previous work. This allows us to reconstruct changes in firn thickness and temperature difference at the South Pole between 30 and 5 kyr BP. We find that variability in firn thickness is controlled in part by changes in snow accumulation rate, which is, in turn, influenced strongly by the along-flowline topography upstream of the ice core site. Variability in our firn temperature difference record cannot be explained by annual-mean processes. We therefore propose that the ice core gas isotopes contain a seasonal bias due to rectification of seasonal signals in the upper firn. The strength of the rectification also appears to be linked to fluctuations in the upstream topography. As further evidence for the existence of rectification, we present new data from the Dome Fuji ice core that are also consistent with a seasonal bias throughout the Holocene.

Our findings have important implications for the interpretation of ice core gas records. For example, we show that the effects of upstream topography on ice core records can be significant at flank sites like the South Pole – they are responsible for some of the largest signals in our record. Presumably upstream signals impact other flank-flow ice cores such as EDML, Vostok, and EGRIP similarly. Additionally, future work is required to confirm the existence of seasonal rectification in polar firn, to determine how spatially and temporally widespread rectifier effects are, and to incorporate the relevant physics into firn air models.

1 Introduction

Past surface temperatures are commonly inferred from ice cores using the water isotope composition of the ice (δ18Oice), which requires a site-specific calibration of the proxy. Early studies calibrated δ18Oice using its modern-day spatial relationship with mean annual temperature near the ice core site, which was hypothesized to be identical to the relationship with temporal variations in site temperature (e.g., Jouzel et al., 1993). Subsequently, temporal calibrations have become possible for cores from Greenland and Antarctica thanks to the development of independent methods of temperature reconstruction based on borehole thermometry (Cuffey et al., 2016; Dahl-Jensen et al., 1998; Buizert et al., 2021) or gas isotope measurements (Severinghaus and Brook, 1999; Huber et al., 2006; Kindler et al., 2014; Buizert et al., 2014; Orsi et al., 2014; Buizert et al., 2021). Calibrations using these methods have shown that the temporal relationship between gas isotope measurements and temperature can indeed differ significantly from the spatial calibration and can also vary in time. Unfortunately, such independent temperature reconstructions are more challenging for East Antarctic ice cores for two main reasons. First, the low snow accumulation rates at these sites means heat diffusion erases some of the thermal history of surface temperature change that borehole thermometry relies on. Second, the smaller, more gradual surface temperature fluctuations typical of the Antarctic climate result in a lower signal-to-noise ratio for gas isotope thermometry. This has made it more difficult to evaluate the calibration and reliability of δ18Oice as a paleotemperature proxy in East Antarctica.

In this paper, we present the first Antarctic application of gas isotope thermometry with the precision necessary to detect interpretable signals. We describe an improved analytical method for making measurements of the isotopic composition of molecular nitrogen (N2) and argon (Ar) on a single ice core sample and present data from the South Pole ice core between 5 and 30 kyr BP. Our method yields a two-fold improvement in precision compared to previous work, meaning we can measure the isotope ratios with a reproducibility of  0.002 ‰ per atomic mass unit. This allows us to use the isotope measurements to separate the gravitational and thermal components of diffusive fractionation in the firn column and thus quantitatively reconstruct past temporal changes in the height of the diffusive column of firn air and the temperature difference across it. The analytical precision corresponds to an uncertainty of  1 m and 0.3 C in firn thickness and temperature difference respectively. Our measurements span the last glacial period, the deglaciation, and the early Holocene, recording changes in climate and firn properties throughout this time.

This study is important as the most meaningful test yet of gas isotope thermometry in Antarctica. Wider application of comparable, high-precision measurements would provide a benchmark for testing the ability of firn densification models to accurately simulate the thermal properties of the firn column across a wide range of past and present climate conditions and has the potential to improve past temperature reconstructions for the South Pole ice core and at other sites in East Antarctica.

2 Reconstructing firn properties

To reconstruct the firn air diffusive column height (DCH) and vertical temperature difference (ΔTz, the difference between temperature at the surface and the lock-in depth), we measure the isotopic composition of molecular nitrogen and argon (δ15N of N2 and δ40Ar of Ar) in air extracted from ice core samples. All isotope ratios are expressed in delta notation relative to the modern atmosphere in units of per mil (‰).

Importantly, changes in the isotopic composition of atmospheric nitrogen and argon are negligible over the timescales relevant for most ice core studies (< 105 years) (Mariotti, 1983; Sowers et al., 1989; Bender et al., 2008). Therefore, deviations in the ice core gas composition from the modern atmosphere must arise locally in the firn column. Gas transport in the firn is primarily by molecular diffusion, and two processes dominate isotopic fractionation of air: gravitational and thermal fractionation.

In the first case, gravitational settling causes enrichment of the heavy isotopes and molecules at the base of the firn due to the lack of turbulent mixing of the air. The amount of enrichment is described by the following barometric equation (Craig et al., 1988; Sowers et al., 1989; Schwander et al., 1993).

(1) δ grav = exp Δ m a / b g z R T - 1 1000 Δ m a / b g z R T 1000

Here, δgrav is the isotopic deviation (in units of ‰), Δm is the mass difference between the isotope pair a and b, g is the gravitational acceleration, z is the firn air diffusive column height, R is the ideal gas constant, and T is the average temperature of the firn column in Kelvin. It is often useful to make a linear approximation to the exponential (via the first order Taylor expansion), as shown in Eq. (1), which adds a relative error of less than 0.5 % for the range of values considered here. Firn thickness depends on the balance between the rates of snow accumulation and densification with both low temperatures and high accumulation rates resulting in a large z. Because site temperature and accumulation rate are strongly and positively correlated in the climate system, variations in z tend to be muted. Broadly, the spatial pattern across Antarctica is one of thicker firn columns in colder locations, suggesting a dominance of the temperature effect. However, in comparing last glacial maximum (LGM) and pre-industrial values of z in central Antarctica, we find thinner firn columns during the colder LGM (Landais et al., 2006), suggesting a dominance of the accumulation effect (Buizert, 2021).

Second, gas composition is fractionated by temperature gradients in the firn, with heavier isotopes and molecules concentrated at the cold end of the gradient by thermal diffusion fractionation (Severinghaus et al., 1998). The magnitude of the fractionation is given by

(2) δ therm = Ω a / b Δ T z ,

where Ω is the empirically measured thermal diffusion sensitivity of the isotope pair a and b, and ΔTz is the temperature difference between the top and bottom of the diffusive column of air. Positive values of ΔTz and δtherm correspond to the top of the firn column being warmer than the base. At the South Pole, the vertical temperature profile depends broadly on the balance between the downward advection of cold ice from the surface and the upward conduction of geothermal heat. Perturbations to either the mean annual surface temperature, basal geothermal heat flux, ice thickness, or vertical velocities can all influence the firn temperature gradient. The height of the firn column at the South Pole ( 120 m) makes it particularly well suited for recording thermal perturbations because the thermal relaxation time of the firn column scales with the square of the firn column height (Cuffey and Paterson, 2010).

By measuring two isotope ratios, we can mathematically solve for the two components of fractionation, allowing us to calculate the height of the past diffusive column of air and the vertical temperature difference across it. To do so, we first use Eqs. (1) and (2) to express δ15N and δ40Ar as the sum of their respective gravitational and thermal components.


Severinghaus et al. (2003) take advantage of the fact that (in the linear approximation) the gravitational fractionation term is 4 times larger for δ40Ar than for δ15N and define δ15Nexcess, a second-order isotope parameter proportional to the temperature difference.

(5) δ 15 N excess = δ 15 N - 1 4 δ 40 Ar = Ω 15 / 14 - 1 4 Ω 40 / 36 Δ T z

Similarly, we can define δ40ArDCH, a second-order isotope parameter directly proportional to the diffusive column height.

(6) δ 40 Ar DCH = δ 40 Ar - Ω 40 / 36 Ω 15 / 14 δ 15 N = g z R T 4 - Ω 40 / 36 Ω 15 / 14

The final step is to convert from the isotope parameters to ΔTz, the firn temperature difference, and z, the diffusive column height by rearranging Eqs. (5) and (6).


This conversion from isotope ratios to the firn physical properties assumes that the isotope ratios occluded in bubbles at the base of the firn column are in diffusive equilibrium with the local environment and that the only fractionating processes occurring are gravity and thermal gradients. This is generally true for the firn column at an ice core site, although we discuss in Sect. 5.2.4 reasons why this might not be the case at the South Pole, Dome Fuji, and potentially other ice core sites.

3 Methods

3.1 Sample recovery and storage

The South Pole ice core SPC14 (hereafter SPICEcore) was drilled between 2014 and 2016 at a site close to the Amundsen–Scott South Pole Station (Casey et al., 2014; Winski et al., 2019; Epifanio et al., 2020; Souney et al., 2021). Ice cores were transported to the National Science Foundation's Ice Core Facility (NSF-ICF) in Denver, Colorado, where 200 g samples were cut and shipped to Scripps Institution of Oceanography in La Jolla, California. The ice was kept colder than 25 C from coring to analysis.

3.2 Gas extraction and mass spectrometry

Our method for the extraction and purification of the trapped gases is similar to that described by Kobashi et al. (2008) and Orsi (2013). Briefly, an 80 g piece of ice is melted in an evacuated vessel, and the gases are stirred out of solution by a magnetic stir bar. Oxygen is removed by reaction with copper turnings heated to 500 C to prevent interference of the 18O18O isotopologue with the 36Ar beam and to improve 29N2 beam stability. Other interfering gases, such as water vapor and carbon dioxide, are removed by a series of glass u-traps immersed in liquid nitrogen at 77 K, and the remaining nitrogen and argon are cryogenically trapped in a stainless-steel tube immersed in liquid helium at 4 K. The dip tube is removed from the liquid helium and allowed to thaw and re-equilibrate for a minimum of 12 h before being analyzed.

Isotopic ratios of nitrogen (29N2/28N2) and argon (40Ar/36Ar), as well as the argon to nitrogen ratio (40Ar/28N2) of the sample gas, are measured on a dual inlet Thermo Finnigan MAT 252 mass spectrometer. Routine laboratory corrections for source pressure imbalance and the Ar/N2 chemical slope are made. Isotope and elemental ratios are expressed in units of per mil (‰) relative to the modern atmosphere, sampled in La Jolla, California, USA. The La Jolla air samples are processed similarly to ice samples, meaning that small biases induced by gas handling cancel out to first order.

We make two important modifications to the methods described by Kobashi et al. (2008) and Orsi (2013). The first is a chemical slope correction to δ15N, which is artifactually enriched by the presence of H2 in the sample gas. The second is the inclusion of a 30 min delay between admission of the sample and reference gas into the bellows and the beginning of the measurement sequence. This is necessary due to an initial measurement bias caused by cooling of the bellows during expansion of the reference gas, which is at a higher pressure than the sample gas prior to expansion. Both modifications are discussed in more detail in Sect. S1 in the Supplement.

3.3 Firn densification modeling

In this work we perform firn densification model simulations using a coupled firn densification and heat transport model that has been described previously elsewhere (Buizert et al., 2014, 2021). The model uses Herron–Langway densification physics formulated in terms of overburden pressure to allow for non-steady-state conditions (Eq. 4c in Herron and Langway, 1980). Firn thermal conductivity is based on Calonne et al. (2019), and other firn and ice thermal properties are based on Cuffey and Paterson (2010). The forward model is forced using the surface temperature and accumulation rate histories at the site. The model simulates the time evolution of firn density and temperature with depth. The close-off density is estimated using the parameterization of Martinerie et al. (1994). Ice core gas properties (gravitational and thermal fractionation and gas age–ice age difference) are calculated and saved at the lock-in density, which is determined using the established approach by Blunier and Schwander (2000) of finding the lock-in density by subtracting a constant value from the Martinerie close-off density. Blunier and Schwander recommend a constant value of 14 kg m−3 at Summit, Greenland. We use a value of 15 kg m−3 based on modern-day observations at the South Pole. The larger value reflects the fact that the South Pole has a very thick lock-in zone. The DCH is equal to the lock-in depth minus the convective zone thickness; thermal fractionation is calculated using the temperature difference between the bottom of the convective zone and the lock-in depth. The convective zone thickness is set to 6 m and the firn surface density to 380 kg m−3 following observations (Todd A. Sowers and Christo Buizert, personal communication, 2021). Ice thickness and geothermal heat flux are held constant at 2600 m and 56 W m−2 respectively – these values are not very well known as SPICEcore was not drilled to bedrock. The model can be run in an inverse mode, in which an automated algorithm is used to estimate the temperature and accumulation histories that best fit the observational δ15N data and the empirically reconstructed estimates of the gas age–ice age difference (Epifanio et al., 2020). We will refer to the optimized inverse scenario as the reference (REF) run; we later describe various model experiments that deviate from the REF scenario.

4 Results

We analyzed samples from 170 depths in SPICEcore between 490 and 1310 m depth. The samples encompass bubble ice, clathrate ice, and the transition zone, where bubbles and clathrates coexist. We measured 14 depths in duplicate, giving us an estimate of analytical reproducibility. Our samples cover the time period from approximately 5000 to 30 000 yr BP at an average resolution of 150 years on the SP19 gas chronology (Epifanio et al., 2020). The measurements were made in two periods, between January and April 2018 and between October and December 2018. We calculate gravitationally corrected δAr/N2 (δAr/N2 grav) by making the common assumption that the enrichment per mass unit is equal to the measured δ15N value (Craig et al., 1988; Bender et al., 1995) (thermal fractionation is negligible compared to the precision of the δAr/N2 measurement). We also make a small gas loss correction to δ40Ar based on δAr/N2 grav, the details of which are described in Sect. S2.

4.1 Reproducibility

We assess the reproducibility of our data by calculating the pooled standard deviation, spooled, which allows us to combine our replicate measurements and evaluate their deviations from their respective means:

(9) s pooled = j = 1 m i = 1 n j δ i , j - δ j 2 j = 1 m n j - m 1 2 ,

where δi,j is the ith delta value for a replicate sample from the jth depth, δj is the mean for all replicate samples for a given depth, nj is the number of samples analyzed for a given depth, and m is the number of depth means analyzed.

Five separate flasks of La Jolla air were analyzed between 5 and 11 times each with at least one flask measured at the start and end of each measurement period. The pooled standard deviation of δ15N, δ40Ar, δAr/N2 grav, and δ15Nexcess for the 40 total La Jolla air measurements is shown in Table 1. We achieve a two- and three-fold improvement, relative to Kobashi et al. (2008) and Orsi (2013), for ice measurements of δ15N and δ15Nexcess respectively. We also note smaller improvements in the reproducibility of the other measurements and that some of the improvement may be due to superior ice quality for SPICEcore. This advance in measurement precision makes it possible to reliably detect the δ15Nexcess record of climatic signals in Antarctic ice for the first time.

Table 1Mass-normalized pooled standard deviation of replicate measurements of δ15N, δ40Ar, δAr/N2 grav, and δ15Nexcess from either reference gas runs (REF), La Jolla air flask samples (LJA), South Pole ice core samples (SPC), or other ice core samples. Units for all four isotope ratios are per mil per atomic mass unit (‰ amu−1), and the mass differences are 1, 4, 12, and 1 amu respectively. The final column indicates n, the number of samples used in the calculation.

Download Print Version | Download XLSX

It is also noteworthy that the mass-normalized pooled standard deviation of δ15Nexcess is smaller than that of δ15N and δ40Ar for the La Jolla air flask (LJA) and South Pole ice core (SPC) samples. This suggests that the measured isotope ratios contain some mass-dependent variability that cancels out when we calculate δ15Nexcess. The reproducibility of the reference gas samples does not show the same pattern, suggesting that the variability is introduced to the LJA samples during gas extraction rather than the mass spectrometry. For the SPC samples, another possibility is that the pattern is caused by real mass-dependent variability in the ice due to well-documented spatial heterogeneity in the depth of bubble close-off on a horizontal length scale of a few centimeters, i.e., similar to the width of an ice core sample (Orsi, 2013). This highlights the importance of measuring δ15N and δ40Ar on the same piece of ice. If δ15N and δ40Ar were measured on different pieces of ice, even adjacent pieces from the same depth in the core, this variability would not cancel out and would increase the scatter in δ15Nexcess.

Finally, we note that the pooled standard deviation of δAr/N2 grav is much worse for the ice samples compared to the LJA measurements. This is because of similar centimeter-scale spatial heterogeneity in argon gas loss during bubble close-off and sample storage. Adjacent pieces of ice are likely to have lost different amounts of Ar and so would not be expected to have the same δAr/N2 grav value.

4.2 Isotope and elemental ratios (δ15N, δ40Ar, and δAr/N2)

Our isotope ratio measurements are shown in Fig. 1a. There is a strong positive correlation (r=0.991) between δ15N and δ40Ar/4 as the variability in both is dominated by gravitational fractionation, which affects mass-normalized isotope ratios equally. As described above, the difference between δ15N and δ40Ar/4, termed δ15Nexcess, reflects thermal fractionation in the firn column. δ15N ranges from a minimum of 0.492 ‰ to a maximum of 0.626 ‰ with a mean of 0.562 ‰. δ40Ar/4 has a range of 0.497 ‰ to 0.625 ‰ with a mean of 0.569 ‰. Temporal variations are discussed below in the context of the firn properties calculated from the isotope ratios.

Gravitationally corrected δAr/N2 (δAr/N2 grav) is depleted relative to the modern atmosphere over much of the depth range of our measurements, with values as low as 6.3 ‰. This is typical of ice core gas records and is due to preferential loss of Ar from the ice during bubble close-off (Craig et al., 1988; Bender, 2002; Severinghaus and Battle, 2006). However, there is also an interval of elevated δAr/N2 grav values between 8 and 18 kyr BP, with values as high as 7.8 ‰. This corresponds to the bubble–clathrate transition zone (BCTZ), where gas molecules are held in coexisting bubbles and clathrate hydrates. Here, δAr/N2 grav is enriched because post-coring gas loss occurs primarily from the bubbles, which are enriched in N2 due to the stronger affinity of Ar for the clathrate phase (Bender et al., 1995; Ikeda-Fukazawa et al., 2001). The transformation to clathrates occurs heterogeneously throughout the core, increasing the scatter in our δAr/N2 grav measurements. Both the enrichment and increased scatter of elemental ratios in the BCTZ have been noted in many ice cores (Suwa and Bender, 2008a, b; Kobashi et al., 2008; Shackleton et al., 2019), but recent work appears to confirm that there is no appreciable isotope fractionation associated with clathration (Oyabu et al., 2021).

Figure 1(a) SPICEcore measurements of δ15N and δ40Ar. The δ40Ar data are divided by 4 so that they can be plotted on the same axis as δ15N, with the visual offset between the two isotope ratios equal to δ15Nexcess (or, equivalently, ΔTz). (b) The firn diffusive column height (DCH) and (c) temperature difference (ΔTz) and equivalent isotope parameters (see Sect. 2 for explanation). (d) SPICEcore measurements of δAr/N2 after correction for gravitational fractionation. All data are plotted on the bottom x axis on the SPC19 gas chronology. The corresponding depths in SPICEcore are indicated on the top x axis. Error bars in DCH and ΔTz represent 1 pooled standard deviation of replicate samples. Error bars for δ15N, δ40Ar/4, and δAr/N2 grav are smaller than the data markers.


4.3 Reconstructed firn properties

We calculate DCH and ΔTz from our isotope data using Eqs. (7) and (8). The time series are shown in Fig. 1. Both DCH and ΔTz increase over the course of the record with DCH increasing from a glacial average of 103 to 111 m in the Holocene and ΔTz increasing from 1.9 to 0 C.

The minimum DCH we reconstruct is  95 m and occurs around 23 kyr BP. The maximum DCH is  115 m at  11 kyr BP, although DCH is also > 110 m around 20 kyr BP and for much of the deglaciation and early Holocene. The minimum (i.e., most negative) ΔTz we reconstruct is 4 C, which occurs around 20 kyr BP, concurrently with a local maximum in firn thickness. In fact, this inverse relationship is a persistent pattern on timescales of a few millennia. Despite a broad positive correlation between DCH and ΔTz over the full 25 kyr record, a negative relationship exists between the higher-frequency variability for much of the record (Fig. 2). This is most evident prior to  17 kyr BP, when the fluctuations in DCH and ΔTz are the largest in amplitude and are clearly inverse of one another but also exist in the younger part of the record (< 12.5 kyr BP). For example, the most positive values of ΔTz around 7 kyr BP are associated with a local minimum in DCH. The slope of the relationship is similar for glacial and Holocene samples, implying that the same physical process may be responsible. During the deglaciation, the inverse relationship between DCH and ΔTz breaks down, with both properties increasing through time. This time period is responsible for the overall positive correlation between the two time series.

Figure 2Firn temperature difference plotted against the diffusive column height. Error bars represent 1 pooled standard deviation of replicate samples. In panel (a), the data are plotted in grey along with a regression to the entire dataset. In panel (b), data older than 17 kyr BP are plotted in blue, data younger than 12.5 kyr BP are plotted in red, and a least squares linear fit to each subset is shown. The slope and squared correlation coefficient, r2, of each fit is also indicated.


5 Discussion

The overall increase in both DCH and ΔTz through the deglaciation is the expected response to increased snowfall due to a warming climate. DCH depends strongly on accumulation rate, which is higher in warmer interglacial periods than in glacials. For ΔTz, low accumulation in the glacial period results in negative values as geothermal heat warms the base of the firn. Then, as the accumulation rate increases, greater downward advection of cold surface ice makes the firn column closer to being isothermal, making ΔTz less negative. The concurrent increase in surface temperature during the deglaciation itself also acts to make ΔTz more positive. Next, we discuss in more detail the processes responsible for the higher-frequency variability in DCH and ΔTz.

5.1 DCH (diffusive column height)

First, we consider the mechanisms that drive changes in DCH. Winski et al. (2019) previously presented a record of SPICEcore δ15N, which, in the absence of complementary δ40Ar measurements, they interpret solely as a gravitational fractionation/firn thickness signal. They argue, based on firn modeling experiments, that firn thickness variability in the Holocene is primarily controlled by the local accumulation rate. Our work supports this interpretation where we have data (11–5 kyr BP). Most of the variability in our δ15N data is due to changes in gravitational fractionation, and we note the correspondence between DCH and the SPICEcore record of accumulation between 11 and 5 kyr BP (Fig. 3).

To further develop our understanding of the mechanisms driving changes in DCH, we must therefore consider the mechanisms that drive changes in the accumulation rate. In SPICEcore, Holocene accumulation rate variability is almost entirely explained by the spatial variability in accumulation upstream from the SPICEcore site (Lilien et al., 2018). This is because the South Pole is located far from an ice divide, with ice flowing at 10 m yr−1 in the direction of 40 W (Hamilton, 2004; Casey et al., 2014). Therefore, the snow deposition site for SPICEcore ice is further upstream for deeper, older ice. In this way, spatial accumulation variability is recorded as temporal variability in the ice core as more distant spatial anomalies are advected to a greater depth below the present-day SPICEcore site. The upstream spatial variability is in turn controlled directly by the local topography (Hamilton, 2004; Fudge et al., 2020). Namely, there is a positive correlation between the accumulation rate and the topographic curvature (second derivative along the direction of the flowline) (Fig. 3, r=0.55). The relationship is evident for at least 100 km in the upstream direction from the South Pole and is consistent with findings at other sites in Antarctica (Waddington et al., 2007; Leonard et al., 2004) and Greenland (Miège et al., 2013; Hawley et al., 2014). The mechanism is that katabatic winds accelerate down slopes as the topography becomes steeper and decelerate as it becomes less steep. This results in greater erosion of snow from ridges (negative second derivative of elevation) and greater deposition in depressions (positive second derivative of elevation). In sum then, Holocene DCH is controlled in part by the upstream topography via its dependence on the accumulation rate. This is most evident in our data between 8.5 and 6.5 kyr BP, when an  8 m local minimum in DCH is co-located with a minimum in the modern spatial pattern of accumulation and with the steepest topographic slope upstream of SPICEcore.

The comparison between spatial (upstream) and temporal (SPICEcore) variability is less straightforward prior to 10 kyr BP because the exact position of the flowline is less certain and changes in climate are expected (Fudge et al., 2020). However, we hypothesize that the Holocene pattern also operated during the glacial period. For example, between 90 and 100 km upstream of the South Pole, the topographic slope is close to or less than zero for the only extended period in the survey data (Fig. 3). The survey line from Lilien et al. (2018) terminates at 100 km, but data from the PolarGAP airborne radar campaign (Jordan et al., 2018) confirm that this feature is part of a broad topographic low on the flank of Titan Dome. The low is associated with a local maximum in both topographic curvature and upstream accumulation, suggesting that the topography does indeed cause higher accumulation in this region in modern times. Using a likely flowline, we find that ice of 20 kyr BP age would have originated at this topographic low. We therefore argue that the local maximum in DCH at 20 kyr BP is due to greater net accumulation in the topographic low (Fig. 3). Other Antarctic ice core records of δ15N and accumulation rate show no evidence for continent-wide climatic changes at this time (Buizert et al., 2021), supporting our argument that this is a local signal, not a climatic one. Our argument requires certain features of the present-day topography to be unchanged over the past 25 kyr. This is certainly possible if these features are linked to the bedrock topography, as has been documented elsewhere in Antarctica (De Rydt et al., 2013).

Whilst some of the variability in DCH almost certainly reflects climatic changes associated with the deglaciation, it is not surprising that the effects of upstream variability are also present given the location of the SPICEcore site far from a dome. Our work shows that the signals associated with upstream effects can be substantial – the feature between 23 and 18.5 kyr BP is the largest in our record – and emphasizes that caution must be applied when interpreting temporal changes in firn thickness in SPICEcore or other cores from flank sites such as EDML, Vostok, and EGRIP.

Figure 3From top to bottom: (a) surface elevation profile (meters above sea level) along the flowline upstream of the SPICEcore site. The light blue line corresponds to snowmobile-mounted GPS data from Lilien et al. (2018). The dark blue line corresponds to PolarGAP airborne radar data from Jordan et al. (2018). The gridded radar data were interpolated to the flowline using inverse-distance-weighted interpolation. (b) Surface slope (first derivative) and (c) curvature (second derivative) of the elevation profiles, calculated in the direction of the flowline. Colors are as in (a). Curvature of the PolarGAP data is not shown due to the coarser spatial resolution of this dataset. (d) The pink line shows the accumulation rate along the flowline upstream of the SPICEcore site from Lilien et al. (2018). The purple line shows the SPICEcore accumulation rate derived from strain-corrected annual layer thicknesses (Winski et al., 2019). (e) Firn diffusive column height and (f) temperature difference calculated from our isotope data. Upstream data in panels (a)(d) are plotted on the bottom x axis as functions of the transit time to the SPICEcore site based on the “mean” scenario in Fudge et al. (2020). The corresponding distances are shown on the top x axis. SPICEcore data in panels (d)(f) are plotted as functions of age on the SPC19 chronology. Grey shading highlights the changes in DCH and ΔTz between 23 and 18.5 kyr BP and between 8.5 and 6.5 kyr BP, which are discussed in the text.


5.2ΔTz (top-minus-bottom firn temperature difference)

The variability in our record of ΔTz is initially challenging to explain. We would have anticipated a positive correlation between DCH and ΔTz since an increase in the accumulation rate ought to result in a thicker firn column and a weaker influence of geothermal heat on the temperature at the lock-in depth. However, although DCH and ΔTz both increase over the course of the deglaciation, we instead observe a negative correlation between DCH and ΔTz throughout most of the record (Fig. 2). There must be some other mechanism that links variability in ΔTz to either changes in accumulation or the local topography.

Furthermore, the 4 C difference in temperature between the top and bottom of the firn column is much larger in magnitude than the present-day temperature difference at the South Pole, which is approximately 0 C (accumulation = 8 cm yr−1). Other sites on the East Antarctic plateau with present-day accumulation rates comparable to the estimated glacial value at the South Pole (4 cm yr−1) also have smaller values of ΔTz, for example, 0.8 C at Dome C (2.5 cm yr−1) and 1.3 C at Dome Fuji (2.5 cm yr−1). Details of how present-day values of ΔTz were determined are described in Sect. S3. Below, we examine several mechanisms that could explain the extreme negative values of ΔTz and its inverse relationship with DCH in SPICEcore.

5.2.1 Surface temperature change

First, we investigate whether the SPICEcore ΔTz reconstruction can be explained by variations in surface temperature. Changes in mean annual site temperature affect the firn temperature difference as the surface snow warms or cools and the vertical temperature profile in the ice sheet adjusts to a new equilibrium. Surface temperature change might also explain the negative relationship between DCH and ΔTz in our data, with a surface cooling trend typically resulting in a more negative ΔTz and a thicker firn column via reduced densification rates (Herron and Langway, 1980). We estimate the surface temperature history using the dynamical firn densification–heat transport model (Sect. 3.3) in an inverse mode. Briefly, the model adjusts initial estimates of past surface temperature and accumulation rate to best fit proxy-based reconstructions of firn thickness and gas age–ice age difference (Δage). This REF model run is able to produce a good fit to the proxy-based estimates of firn thickness and Δage for SPICEcore (Buizert et al., 2021) and also agrees well with our estimates of DCH and ΔTz for much of the last glacial and Holocene periods (Fig. 4). However, the modeled ΔTz from the REF run does not agree well with our ΔTz reconstruction during the LGM and for much of the deglaciation. The model does not reproduce the most negative values of ΔTz between 23 and 18.5 kyr BP, nor does it capture some of the most positive values between 8.5 and 6.5 kyr BP. The firn model is not capable of fitting the ΔTz data while simultaneously fitting the observational DCH and Δage data.

To evaluate what temperature history would be required to fit the ΔTz data, we perform an additional experiment in which the firn temperature is decoupled from the firn densification physics (we call this the DECOUPLE run). In the DECOUPLE run the firn densification rates are not calculated, but rather they are read out from a data file corresponding to the densification rates from the REF experiments. Accumulation rates are likewise equal to those from the REF run. The inverse model is then tasked to reconstruct the surface temperature history that best fits the ΔTz data. The DECOUPLE surface temperature history required to fit our ΔTz data is show in Fig. 4. Note that the DECOUPLE scenario is not internally consistent as the firn densification rates are inconsistent with the temperature forcing used. The design of the DECOUPLE experiment simply allows us to control the thermal gradient in the firn column (and thus ΔTz) while simultaneously ensuring we use the correct firn thickness and rate of downward advection of ice and heat.

We focus our interpretation of the DECOUPLE simulation on the direction and timing of changes in the inferred temperature history, rather than the absolute values, as ΔTz is more sensitive to changes in surface temperature than to the temperature itself. Also shown for comparison is the optimal temperature history from the REF run (Buizert et al., 2021) and a temperature history from Kahle et al. (2021) based on a calibration of δ18Oice using the SPICEcore Δage data and the diffusion length of water isotopes in the firn. Note that both temperature histories are partially constrained by the Δage data, so they are not wholly independent.

The DECOUPLE temperature history differs substantially from the other temperature estimates. It features a prolonged 5 C cooling between 23 and 20 kyr BP, followed by a rapid 5 C warming from 13 to 11 kyr BP. The cooling event in the decoupled temperature history happens at a time when the other estimates indicate either constant temperatures or a slight warming associated with the initiation of the deglaciation, whereas the decoupled history shows almost no warming until the deglacial temperature change is almost fully realized in the other estimates. The timing and sign of temperature changes in the decoupled temperature history also bear little resemblance to other Antarctic ice core temperature reconstructions (Buizert et al., 2021; Uemura et al., 2018; Cuffey et al., 2016; Stenni et al., 2011; Jouzel et al., 2007; Petit et al., 1999).

Based on both the poor agreement of the measured and modeled ΔTz in the REF run and the poor agreement of the DECOUPLE temperature history with other reconstructions, we conclude that surface temperature change is unlikely to fully explain our record of ΔTz, particularly the most negative values at 20 kyr BP.

Figure 4(a) Modeled Δage at SPICEcore in the REF (light grey) and DECOUPLE (dark grey) experiments. Markers show empirical Δage data from Epifanio et al. (2020). (b) SPICEcore ΔTz reconstruction (green points), together with the modeled ΔTz from the REF (light grey line) and DECOUPLE (dark grey line) model experiments. (c) South Pole surface temperature reconstructions from the REF (light red) and DECOUPLE (dark red) model experiments. Shading shows 1σ model uncertainty from Buizert et al. (2021). Also shown is a temperature reconstruction from Kahle et al. (2021) (orange), based on a calibration of δ18Oice using Δage and the water isotope diffusion length proxy. Shading shows 1σ standard deviation of model ensemble runs.


5.2.2 Ice thickness

Second, we investigate whether the SPICEcore ΔTz reconstruction may be explained through variations in the thickness of the ice sheet. Ice thickness influences ΔTz by controlling the vertical strain rate in the ice sheet and thereby the downward advection of cold ice and the ability of geothermal heat to warm the base of the firn column. Temporal changes in ice thickness over the course of our record are certainly plausible, especially given that the flank location of the South Pole means older ice originated upstream from the present-day SPICEcore site. The relevant parameter for our record of ΔTz is the thickness of the ice column at the time and location that the bubbles were occluded. Variations in ice thickness experienced by SPICEcore ice are therefore the result of both temporal fluctuations in ice sheet elevation and upstream spatial fluctuations in ice thickness.

Temporal changes in ice sheet elevation at the South Pole were estimated by Fudge et al. (2020) using output from a full ice sheet model from Pollard et al. (2016). They conclude that changes in ice thickness have most likely been smaller than ±100 m in the past 20 kyr, with a mean elevation change of +16 m at 20 kyr BP. Certain combinations of model parameter values give changes between 325 and +250 m. Outside of the Holocene, the rate of change is always less than 20 m kyr−1.

Spatial changes in ice thickness upstream of the SPICEcore site are dominated by changes in bed topography. For the Holocene, changes are well-known as the flowline is tightly constrained for this time period and the bed topography has been determined by ice penetrating radar (Lilien et al., 2018). Fluctuations are smaller than ±250 m. Beyond the Holocene, estimates of the bed topography do exist, but the exact flowline position becomes increasingly uncertain. Kahle et al. (2021) offer some constraints from their estimate of the SPICEcore thinning function and infer possible changes in bed elevation of around +200 m between 33 and 26 kyr BP and 200 m between 23 and 18.5 kyr BP.

We therefore seek to determine the maximum plausible change in ΔTz by simulating a 500 m thickening or thinning using a one-dimensional ice flow model (Fudge et al., 2019). A 500 m change is 80 m larger than the range of present-day ice thickness fluctuations in the 100 km upstream of the SPICEcore site, corresponding to the past 20–25 kyr (Lilien et al., 2018); it is more than twice the magnitude of the changes in bed topography inferred by Kahle et al. (2021) in the past 30 kyr BP, and it is 1.5 to 2 times larger than the most extreme modeled surface elevation changes in the past 20 kyr. In our experiments, the modeled changes in ice thickness occur linearly over 2000 years, and the post-change thickness is set to the present-day value of 2800 m in each case. We repeat each experiment twice with accumulation rates of 2 and 4 cm yr−1.

The results of these simulations are shown in Fig. 5. The thinning (thickening) experiments result in a decrease (increase) in ΔTz to more (less) negative values. Lower accumulation rates correspond to more negative values of ΔTz due to weaker downward advection of surface heat. The change in ΔTz is less than 1 C for each experiment – 3 times smaller than the largest change in our record. We therefore conclude that fluctuations in ice thickness are also unable to fully explain our observations.

Figure 5Modeled firn temperature difference for a thickening (blue lines) or thinning (red lines) of the ice column. Solid lines correspond to experiments with an accumulation rate of 2 cm yr−1, and dashed lines correspond to an accumulation rate of 4 cm yr−1.


5.2.3 Basal geothermal heat flux

Third, we investigate whether the SPICEcore ΔTz reconstruction may be explained through variations in the basal geothermal heat flux (GHF). As discussed above, most East Antarctic ice cores have negative values of ΔTz caused by geothermal heat impinging on the base of the firn column due to low accumulation rates at these sites. Although the GHF at the South Pole is constrained by borehole temperature measurements (Price et al., 2002; Beem et al., 2018), the firn temperature gradient at 20 kyr BP may have been set by very different basal conditions as the ice sheet flowed over regions of greater or lesser GHF towards the present-day SPICEcore site. In fact, a recent survey upstream of the SPICEcore site inferred values as high as 120 W m−2 due to local faulting and hydrothermal activity (Jordan et al., 2018), more than double previous estimates for the region from continent-scale models (Van Liefferinge and Pattyn, 2013).

To test the hypothesis that the most negative values of ΔTz and the negative relationship between DCH and ΔTz are the result of spatiotemporal variations in GHF, we simulate the effect of a step change in the GHF using the firn densification model described in Sect. 3.3. To calculate an upper bound on the plausible change in ΔTz, we choose a low starting value of 40 W m−2 and either double or triple the GHF instantaneously. We repeat the experiments with three difference ice thicknesses: 2800, 2300, and 1500 m. The first represents the present-day ice thickness at the South Pole, and the second represents a plausibly thinner ice sheet (as deduced in Sect. 5.2.2). The third is a more extreme scenario, representing the minimum observed thickness in a recent survey of the upstream region (Beem et al., 2021). For all experiments, the accumulation rate is 4 cm yr−1, and the surface temperature is 58 C, representing LGM conditions.

The results of the model experiments are shown in Fig. 6. Over the course of  100 kyr, the firn temperature profile adjusts to the new steady state, with ΔTz decreasing by 0.7–2.9 C, depending on the ice thickness and GHF change. However, in our record of ΔTz, variability of this magnitude occurs in  5 kyr. The model indicates that GHF changes can explain < 0.5 C of change in ΔTz on this timescale and only under the most extreme case of a tripling of GHF with a 1500 m ice column. Furthermore, a larger GHF would likely result in a slightly smaller DCH as the warmer firn column would densify more rapidly. This prediction is confirmed by the model (not shown) and is opposite to the inverse relationship we observe between DCH and ΔTz. Therefore, we conclude that spatial variability in the GHF upstream of the South Pole is unable to fully explain our observations, particularly the rapid 3 C changes in ΔTz between 23 and 18.5 kyr BP.

Figure 6Modeled change in ΔTz in response to a change in the basal geothermal heat flux (GHF), plotted as the difference from the value at t=0. Solid lines correspond to experiments in which the GHF is doubled from 40 to 80 W m−2, and dashed lines correspond to experiments in which the GHF is tripled from 40 to 120 W m−2. Each experiment is repeated for three different ice thicknesses: 1500, 2300, and 2800 m. The vertical black line indicates t=0, the time at which the change in GHF occurs in the model. The box in (a) indicates the area covered by (b).


5.2.4 Rectification of seasonal thermal signals

Last, we investigate whether the SPICEcore ΔTz reconstruction may be explained through so-called rectifier effects. Based on the evidence presented in previous sections, we argue that none of the processes known to control the annual-mean firn temperature profile can adequately explain our ΔTz observations in terms of their magnitude, rate of change, and inverse relationship with DCH. Notably, however, temperature differences much larger than 4 C do arise on sub-annual timescales within the upper 20 m of the firn, at the South Pole and elsewhere, in response to the seasonal surface temperature cycle (Dalrymple et al., 1966; Severinghaus et al., 2001; Brandt and Warren, 1997; Town et al., 2008). The corresponding gas isotope thermal fractionation signals only penetrate to  10–15 m depth (Severinghaus et al., 2001; Weiler et al., 2009) and are typically assumed to cancel out each year such that the deep firn (and thus the ice core gas archive) reflects the annual mean. Although existing firn air data from multiple sites are largely consistent with this assumption, the data are often lower precision than our measurements and are unlikely to represent a complete picture of firn processes on the spatial and temporal scales captured by our ice core record. Therefore, because annual-mean processes are unable to explain our data, we now investigate the possibility that some values of ΔTz in our SPICEcore record could be the result of isotope signals in the deep firn being biased towards a particular season at certain times in the past. During winter, cold surface ice overlays a warmer firn column producing negative ΔTz values in the upper firn, and vice versa during summer. To explain the most negative values of ΔTz between 23 and 18.5 kyr BP, we infer a wintertime bias that is either weaker or non-existent at other times in the record. In addition, we propose that the most positive values between 8.5 and 6.5 kyr BP may represent a summertime bias. One notable aspect of this mechanism is that its strength can change very quickly – firn air convection appears and disappears on seasonal timescales. This may help to explain the changes in ΔTz we observe between 23 and 18.5 kyr BP that are either too large or too abrupt to be explained by the processes discussed in previous sections.

In the sections below, we discuss mechanisms that might produce a summer or winter bias and argue that this hypothesis can explain many features of our dataset, including the most negative values of ΔTz, the rate at which they develop, and the inverse relationship between ΔTz and DCH.

Wintertime bias due to Rayleigh–Bénard convection

Seasonal isotopic thermal fractionation signals in the firn are typically overwritten by the opposite signal of the following season (Severinghaus et al., 2001; Weiler et al., 2009). One way a seasonal bias can develop in the deep firn is if one season's isotope signal is preferentially preserved by being advected down into the firn, below the depth to which the next season's diffusive isotope signal penetrates. This type of differential preservation of winter versus summer signals due to covariation of gas transport and concentration has been called a “seasonal rectifier effect” in prior literature (Denning et al., 1995; Severinghaus et al., 2001, 2010; Dreyfus et al., 2010; Trudinger et al., 2020). We adopt this language here.

This type of rectification requires a slow, non-turbulent, downward movement of air that occurs during one season but not the other. A plausible driving mechanism is the snow temperature inversion that arises in winter. Because snow and firn are efficient emitters in the infrared band and are usually warmed from below, their temperature is often coldest at the surface. This is especially true in winter when incoming solar radiation is reduced or even absent. The temperature inversion results in an unstable air density profile in the firn that can trigger buoyancy-driven Rayleigh–Bénard convection, thus advecting seasonal isotope signals deeper into the firn. In this section we discuss evidence for this type of air movement in snow and firn and investigate its ability to explain our SPICEcore gas isotope records.

Sturm and Johnson (1991) demonstrated that buoyancy-driven overturning occurs readily in sub-Arctic snow in Alaska. By making hourly observations of the three-dimensional temperature field within the winter snowpack for 3 years, they were able to observe large horizontal temperature gradients within the snow that were initiated and maintained by columns of rising warm air and sinking cold air. This convection occurred almost continuously throughout two successive winters. There is also ample evidence for air circulation within snow and firn from Antarctica, particularly if vertical cracks allow for fast upward return flow (Giovinetto, 1963; Albert et al., 2004; Fahnestock et al., 2004; Courville et al., 2007; Severinghaus et al., 2010). Unfortunately, direct observations of changes in firn air composition associated with convection are scant since firn air sampling happens almost exclusively in the summer. However, there are published data from a winter firn air sampling campaign at the South Pole. In this case, the authors did indeed find that the peak wintertime isotope signal occurred deeper than their firn air model predicted and speculated that this could be due to downward transport of the isotope anomaly by slowly sinking air (Severinghaus et al., 2001). If correct, this would provide confirmation not only of wintertime convection at the South Pole but also that thermal isotope signals can be carried down into the firn by convection without being destroyed by turbulent mixing.

To test their hypothesis, we compare their wintertime firn air measurements from the South Pole with values predicted by firn air model runs with and without parameterized Rayleigh–Bénard convection (Fig. 7). In the model run without convection, the gases diffuse towards gravitational and thermal equilibrium as they are slowly advected downwards with the densifying firn and occluded in bubbles in the lock-in zone. Because the model is one-dimensional, it is not possible to explicitly simulate a three-dimensional Rayleigh–Bénard convection cell. Instead, we model just the sinking core of a convection cell, which we parameterize as an 8 cm d−1 downward transport of gas between 0 and 20 m. Between 20 and 25 m, the downward transport decays to zero, resulting in mass convergence that would be balanced in the real world by horizontal transport and a return flux of gas to the surface. This approach allows us to approximate how the gas isotopes respond to convection using a one-dimensional model. The model run with downward transport agrees better with the observed wintertime firn air isotope ratios, with the negative wintertime values occurring deeper in the firn than in the model run with no downward advection. The model and the data therefore support our hypothesis that convection can carry seasonal thermal isotope signals down into the firn.

Figure 7Shallow firn air data from the winter sampling campaign at the South Pole by Severinghaus et al. (2001). Solid lines are a firn air model run that includes a slow downward advection of air between 0 and 25 m (see text for details). Dashed lines are the model run without any downward advection, as shown in Fig. 5 in Severinghaus et al. (2001).


Because isotope data are only available in the top 16 m of the firn, we do not have an observational constraint on the strength of rectification in the deep firn where ice core signals are recorded. To demonstrate that seasonal convection can affect isotope values in the deep firn, we perform an additional experiment with the firn air model. We simulate the isotope values in the full firn column under idealized South-Pole-like conditions (110 m thick firn, 51 C annual-mean temperature, 7 cm yr−1 accumulation) and impose a 14 cm d−1 downward advection throughout winter (April–September). In the model, the wintertime signal is advected deeper than the summer signal and so is not fully canceled out. This results in a 0.008 ‰ bias in the annual-mean signal in the deep firn compared to the control run with no downward advection (Fig. 8, Movie S1 in the Supplement). The bias is of comparable magnitude to the signals in our SPICEcore record, demonstrating that this mechanism could plausibly explain some of the millennial variability we observe.

Figure 8Results of idealized modeling experiment. Panels (a) and (b) show the temperature and advection forcing applied to the firn air model. The solid lines correspond to the “with rectifier” run, and the dotted line in (b) corresponds to the “without rectifier” run with no vertical advection. Panels (c) and (d) show the vertical profile of δ15Nexcess in the firn column at the end of summer and winter respectively. The grey line is the run without advection, and the green line is with advection. The days corresponding to the profiles are indicated by the vertical red and blue lines in the upper panels. An animated version of this figure is available in the Supplement as Movie S1.


To explain the correlation between DCH and ΔTz in our SPICEcore record, the strength of the wintertime convection must be linked to the wind speed and/or the topographic slope such that the rectifier is strongest when DCH is also at its maximum. We hypothesize that this link is provided by the energy balance at the snow surface. Stronger katabatic winds on steeper slopes weaken the air temperature inversion by turbulently mixing heat down to the surface from aloft (Hudson and Brandt, 2005; Pietroni et al., 2014). The opposite is true in areas of minimal slope: weaker winds allow a strong inversion to develop via efficient loss of infrared radiation to space from the surface snow. This intense cooling of the surface promotes convection in the firn (Sturm and Johnson, 1991), which would strengthen the wintertime bias. Low wind speeds probably also limit the formation of low-permeability wind crusts that would inhibit convection (Domine et al., 2018). By this mechanism, the wintertime bias would be strongest and ΔTz most negative in areas of flat topography, as we observe in our SPICEcore record (Fig. 3).

As further evidence for this type of seasonal rectifier, we also present a previously unpublished ΔTz record from the Dome Fuji ice core. The core was drilled in 1994–1996, and samples were stored at 50 C until they were analyzed at Scripps Institution of Oceanography in 2007 using a different method to our SPICEcore dataset (Bereiter et al., 2018). Briefly, an ice sample of 800–900 g was melted in an evacuated vessel, and the released air was continuously transferred to a dip tube through a 100 C water trap while stirring the meltwater. The air sample was split in two aliquots (Method 1 in Bereiter et al., 2018): one was measured with Thermo Delta-Plus XP for δ15N, and the other was gettered to extract noble gases and then measured with Thermo Finnigan MAT 252 for δ40Ar. The isotope data and the reconstructed ΔTz data are shown in Fig. 9, and we compare them to our estimate of the modeled Holocene ΔTz from Buizert et al. (2021). The model estimate is based on the same firn densification modeling approach described in Sect. 3.3 and constrained by Dome Fuji δ15N and empirical Δage datasets described in Buizert et al. (2021). To estimate the uncertainty in the modeled ΔTz, we re-run the model with different values of the GHF and accumulation rate. We change the GHF by ±10 W m−2 and the accumulation rate by ±10 %. The total uncertainty we report is the quadrature sum of the difference between these model runs and the optimal scenario.

Just like the SPICEcore record, the Dome Fuji ΔTz data show evidence of a wintertime bias due to rectification. The mean of the Holocene ΔTz data is more negative than both the present-day ΔTz and the modeled Holocene ΔTz. Large changes in surface temperature, ice thickness, and GHF can be excluded during the Holocene, so we conclude that the mismatch is most likely due to rectification producing a wintertime bias throughout the Holocene at Dome Fuji. Because katabatic winds are weak at ice domes due to the flat topography, we expect that the wintertime Rayleigh–Bénard rectifier would be particularly effective at this site. This finding strengthens the case for the existence of rectification in Antarctica and demonstrates that rectification can affect gas records at both dome and flank sites and over a wide range of site characteristics (Dome Fuji is 1000 m higher in elevation, is 5 C colder, and receives half as much snow accumulation).

Also plotted is the ΔTz calculated from δ15N and δ40Ar measurements on firn air collected at Dome Fuji in 1998, which is 1.2 C (Fig. 9, Sect. S3.2). This is more positive than the Holocene ice data and is consistent with the present-day observed firn temperature profile, suggesting no winter rectification is necessary to explain current conditions at Dome Fuji. This could be due to the cessation of rectification at some time during the past 2000 years, perhaps in the last century due to anthropogenic warming (the ice surface absorbs downwelling longwave radiation from greenhouse gases very effectively).

Figure 9Measurements of (a) δ15N and δ40Ar used to calculate (b) δ15Nexcess and an estimate of ΔTz from the Dome Fuji ice core. The ΔTz data are plotted as dark green circles and compared to a model estimate of past ΔTz at Dome Fuji from Buizert et al. (2021) (grey line and shading). The dashed green line shows the mean of the data, and the shading represents 1 standard error of the mean of the six samples. The light green point shows an estimate of modern ΔTz at Dome Fuji calculated using the method described in Sect. S3.2. The estimate is based on a new firn air dataset from archived samples collected in 1998 (Kawamura et al., 2006) and re-measured at SIO in 2008.


In summary, we propose that low wind speeds over areas of minimal topographic slope cause surface snow temperatures to be colder than on steeper slopes. In winter, this can result in an unstable air density profile in the firn and slow, non-turbulent convection of air to a depth of 10–20 m. This is deep enough to produce a cold, wintertime bias in our ice core records of ΔTz. In the Dome Fuji ice core, this bias existed throughout the Holocene until at least 2 kyr BP, whereas in SPICEcore, the cold bias is strongest at 20 kyr BP and is co-located with a thicker firn column due to the increased net accumulation of snow associated with slower and/or decelerating winds. Although this hypothesis is somewhat speculative, we believe this mechanism can plausibly explain (i) the most negative values in our record of ΔTz, (ii) the observed rate of change in ΔTz, and (iii) the inverse relationship with DCH.

Summertime bias due to turbulent convective mixing

The slow, non-turbulent air circulation described above results in a wintertime bias in the deep firn. However, some sites in Antarctica experience vigorous turbulent mixing in the upper few meters of the firn column – termed the convective zone (Sowers et al., 1992; Bender et al., 1994; Kawamura et al., 2006; Severinghaus et al., 2010). This convective mixing of the free atmosphere into the surface firn “resets” the air composition back to atmospheric values, eroding the seasonal, gas isotope thermal fractionation signals. The depth and extent of the mixing is controlled in part by the surface wind speed, with deeper convection associated with faster winds (Kawamura et al., 2006). Because katabatic winds are generally stronger in winter (van den Broeke and van Lipzig, 2003), we propose that a summer bias in ΔTz could originate via a seasonality in the strength of convective mixing in the firn. We might expect stronger wintertime winds to be more effective than summertime winds at eroding the thermal signals in the upper firn, meaning the summertime thermal signal would be preferentially preserved in the deep firn.

Again, we test the plausibility of this hypothesis with the firn air model. For this experiment, we parameterize the convective zone as an eddy diffusivity term in the upper 10 m of the firn that varies seasonally in magnitude – the eddy diffusivity is 5 times larger in winter than in summer. This dampens the wintertime thermal isotope signal, meaning the summer signal is preferentially preserved, resulting in a +0.008 ‰ bias in the annual-mean signal in the deep firn compared to the control run with no seasonal change in the eddy diffusivity (Fig. 10, Movie S2). Again, the bias is of comparable magnitude to the signals in our SPICEcore record, demonstrating that this mechanism could plausibly explain some of the millennial variability we observe.

The summertime bias hypothesis is consistent with our SPICEcore data in that it predicts a deeper, stronger convective zone on the steeper slopes  50 km upstream of the SPICEcore site, where wintertime katabatic wind speeds would be faster (Vihma et al., 2011). This would produce a stronger, more positive bias in this location, potentially explaining the occurrence of positive values of ΔTz (Fig. 3). Although previous authors have speculated that this type of rectification could affect firn and ice core gas records (Severinghaus et al., 2001, 2010; Dreyfus et al., 2010; Petrenko et al., 2013; Verhulst, 2014), observational evidence is limited to one potential site (Law Dome, Antarctica; Trudinger et al., 2020). Future firn air campaigns may help to uncover additional evidence of rectification via seasonal variability in convective strength.

Figure 10As in Fig. 8, panel (a) shows the temperature forcing and panel (b) shows the eddy diffusivity forcing applied to the firn air model. The solid lines correspond to the “with rectifier” run, and the dotted line in (b) corresponds to the “without rectifier” run with no seasonal change in eddy diffusivity. Panels (c) and (d) show the vertical profile of δ15Nexcess in the firn column at the end of summer and winter respectively. The grey line is the run without advection, and the green line is with advection. The days corresponding to the profiles are indicated by the vertical lines in the upper panels. An animated version of this figure is available in the Supplement as Movie S2.


6 Broader implications and future work

Our work demonstrates that gas isotope thermometry can provide meaningful paleoclimate information from Antarctic ice cores. The improved precision of our analytical method allows us to resolve changes in gravitational and thermal fractionation throughout the last deglaciation, although, for SPICEcore, the effects of upstream topography and possible seasonal rectification prevent us from making a surface temperature reconstruction. Application of measurements like ours to other ice cores is necessary to reveal how important these confounding factors are in other ice core gas records.

Rectification of ice core gas records has received limited attention in the literature so far, but our work argues that more careful consideration is necessary. Failure to recognize and account for rectifier effects where they do exist could potentially lead to incorrect temperature estimates.  Fortunately, it is unlikely that rectifier effects would have been significant for previous gas isotope thermometry studies in Greenland (e.g., Kobashi et al., 2007, 2011; Orsi et al., 2014; Landais et al., 2004, 2006; Huber et al., 2006). The presence of rectification via the mechanism we describe likely requires specific surface conditions such as stagnant air and a strong atmospheric temperature inversion. These conditions probably occur rarely on the Antarctic plateau and are even less common in Greenland. To have any effect on the composition of air in the deep firn and closed-off ice they must persist for many weeks or months at a time and reoccur every year for many decades. Furthermore, in the case of Kobashi et al. (2011), agreement between their temperature reconstruction, regional climate model outputs, and modern instrumental records also supports their analysis and interpretation. However, it might be necessary to account for rectifier effects in future gas isotope thermometry studies in Antarctica.

In principle, all gases would be affected by the processes we describe, not just nitrogen and argon. However, it is important to note that, if the seasonal bias we infer in ΔTz is indeed thermal in origin, rectifier effects are likely smaller than typical signals of interest in many common ice core gas proxies. This is because the effect on isotopic and elemental ratios ought to be proportional to the thermal diffusion sensitivities of the gas pair. Thus, for a 3 C bias in ΔTz, the rectification of CO2 concentration, for example, would be less than 0.3 ppmv (Weiler et al., 2009; Leuenberger and Lang, 2002). For δ15N, the bias is approximately 0.014 ‰ for each 1 C of rectification, corresponding to a 3 m bias in the calculated firn thickness. For gases with seasonal cycles in atmospheric abundance, the amount of rectification will be proportional to the amount of mixing or advection and the magnitude of the seasonal cycle rather than thermal diffusion sensitivity. The signal size will therefore be specific to each site but is likely to be only a few percent of the seasonal variability (Trudinger et al., 2020).

In order to interpret ice core gas records accurately, including gas isotope thermometry data, it is crucial to determine the spatial and temporal prevalence of rectifier effects in Antarctica and Greenland and to learn more about the physical processes responsible. Important goals for future work would be to identify clear evidence for contemporary seasonal rectification in deep firn air and shallow ice samples and to determine the link to air transport in the firn and/or local meteorology. The topography upstream from the South Pole would make a promising candidate site. It is possible that rectification will only affect sites with very specific conditions, meaning temperature reconstruction is a simpler task for other Antarctic cores. Alternatively, it may be possible to identify and correct for rectification effects using the isotope ratios of other inert gases such as Ne, Kr, and Xe (e.g., Kawamura et al., 2013). We also show that it is important to consider the effect of changes in basal geothermal heat flux and ice thickness when interpreting gas isotope thermometry data. The magnitudes of these effects are specific to each ice core site and should be considered when choosing candidate cores for gas isotope thermometry.

7 Conclusions

We present a new analytical method for measuring nitrogen and argon isotopes in ice core samples and the first major Antarctic application of gas isotope thermometry with the precision necessary to resolve typical Antarctic climatic signals. We quantitatively separate gravitational and thermal components of isotopic fractionation to reconstruct past changes in the height of, and temperature difference across, the diffusive firn column at the South Pole. We find that both firn thickness and the firn temperature difference are influenced by local topographic variations along the flowline upstream from the ice core site. The impact of topography generates the largest signals in our record, demonstrating that upstream effects must be considered when interpreting similar proxies in ice cores drilled at flank sites. At the South Pole, firn thickness is greater in areas of negligible topographic slope due to greater net accumulation. The firn temperature gradient is also influenced by the topographic slope potentially due to a seasonal rectification caused by the interaction of katabatic winds with surface topography and air in the uppermost firn column. Although we are unable to conclusively determine the origin of the rectifier, we suggest two mechanisms that could plausibly be responsible. Similar evidence for rectification in the Dome Fuji ice core suggests that both dome and flank sites are susceptible. Until now, seasonal rectification has been assumed to have negligible impact on ice core gas records due to limited observational evidence. Our data show that a more careful assessment of rectification is necessary to ensure accurate interpretation of gas isotope thermometry data from Antarctic ice cores, and our new analytical technique can now be deployed to search for this effect at other sites. Determining how widespread rectification is, both spatially and temporally, is crucial if gas isotope thermometry is to be used more widely in East Antarctica.

Data availability

The SPICEcore δ15N, δ40Ar, and δAr/N2 data associated with this study are archived online at the U.S. Antarctic Program Data Center (, Morgan and Severinghaus, 2022) and are available from the corresponding author upon request.


The supplement related to this article is available online at:

Author contributions

All authors contributed to this study. JDM made the SPICEcore ice core gas measurements with input from JPS. CB performed the firn densification modeling. TJF performed the ice flow modeling. KK made the Dome Fuji firn air and ice core gas measurements. JDM, CB, TJF, KK, JPS, and CMT interpreted the results and contributed to the discussion. JDM wrote the paper with contributions from all authors.

Competing interests

The contact author has declared that none of the authors has any competing interests.


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


We thank the U.S. Ice Drilling Program for drilling the ice core; the 109th New York Air National Guard for airlift in Antarctica; NSF's Antarctic Infrastructure and Logistics and Antarctic Support Contractors, the SPICEcore field team, and the members of the South Pole station who facilitated the field operations; and the National Science Foundation Ice Core Facility for ice core processing and archiving. We thank Ross Beaudette for help with the ice analysis, David Lilien, and Ed Brook for their helpful comments and discussion, and Benjamin Hills, Bob Hawley, and Max Stevens for sharing South Pole firn temperature data. For the South Pole gas isotope measurements and firn densification modeling, we gratefully acknowledge funding from NSF (1443710, 1443472, and 1643394). Kenji Kawamura and the Dome Fuji measurements were supported by JSPS and MEXT KAKENHI grants (18749002, 21671001, 15KK0027, and 17H06320). The manuscript was greatly improved thanks to constructive feedback from Markus Leuenberger and one anonymous referee.

Financial support

This research has been supported by the Office of Polar Programs (grant nos. 1443710, 1443472, and 1643394) and the Japan Society for the Promotion of Science (grant nos. 18749002, 21671001, 15KK0027, and 17H06320).

Review statement

This paper was edited by Joel Savarino and reviewed by Markus Leuenberger and one anonymous referee.


Albert, M., Shuman, C., Courville, Z., Bauer, R., Fahnestock, M., and Scambos, T.: Extreme firn metamorphism: impact of decades of vapor transport on near-surface firn at a low-accumulation glazed site on the East Antarctic plateau, Ann. Glaciol., 39, 73–78,, 2004. 

Beem, L. H., Cavitte, M. G. P., Blankenship, D. D., Carter, S. P., Young, D. A., Muldoon, G. R., Jackson, C. S., and Siegert, M. J.: Ice-flow reorganization within the East Antarctic Ice Sheet deep interior, Geol. Soc. Lon. Spec. Publ., 461, 35–47,, 2018. 

Beem, L. H., Young, D. A., Greenbaum, J. S., Blankenship, D. D., Cavitte, M. G. P., Guo, J., and Bo, S.: Aerogeophysical characterization of Titan Dome, East Antarctica, and potential as an ice core target, The Cryosphere, 15, 1719–1730,, 2021. 

Bender, M., Sowers, T. A., and Lipenkov, V.: On the concentrations of O2, N2, and Ar in trapped gases from ice cores, J. Geophys. Res.-Atmos., 100, 18651–18660,, 1995. 

Bender, M. L.: Orbital tuning chronology for the Vostok climate record supported by trapped gas composition, Earth Planet. Sc. Lett., 204, 275–289,, 2002. 

Bender, M. L., Sowers, T., Barnola, J.-M., and Chappellaz, J.: Changes in the O2/N2 ratio of the atmosphere during recent decades reflected in the composition of air in the firn at Vostok Station, Antarctica, 21, 189–192,, 1994. 

Bender, M. L., Barnett, B., Dreyfus, G., Jouzel, J., and Porcelli, D.: The contemporary degassing rate of 40Ar from the solid Earth, P. Natl. Acad. Sci. USA, 105, 8232–8237,, 2008. 

Bereiter, B., Kawamura, K., and Severinghaus, J. P.: New methods for measuring atmospheric heavy noble gas isotope and elemental ratios in ice core samples, Rapid Commun. Mass Spectrom., 32, 801–814,, 2018. 

Blunier, T. and Schwander, J.: Gas enclosure in ice: age difference and fractionation, in: Physics of Ice Core Records, 307–326, 2000. 

Brandt, R. E. and Warren, S. G.: Temperature measurements and heat transfer in near-surface snow at the South Pole, 43, 339–351,, 1997. 

Buizert, C.: The Ice Core Gas Age-Ice Age Difference as a Proxy for Surface Temperature, J. Glaciol., 48, e2021GL094241,, 2021. 

Buizert, C., Gkinis, V., Severinghaus, J. P., He, F., Lecavalier, B. S., Kindler, P., Leuenberger, M., Carlson, A. E., Vinther, B., Masson-Delmotte, V., White, J. W. C., Liu, Z., Otto-Bliesner, B., and Brook, E. J.: Greenland temperature response to climate forcing during the last deglaciation, Geophys. Res. Lett., 345, 1177–1180,, 2014. 

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

Calonne, N., Milliancourt, L., Burr, A., Philip, A., Martin, C. L., Flin, F., and Geindreau, C.: Thermal Conductivity of Snow, Firn, and Porous Ice From 3-D Image-Based Computations, Geophys. Res. Lett., 46, 13079–13089,, 2019. 

Casey, K. A., Fudge, T. J., Neumann, T. A., Steig, E. J., Cavitte, M. G. P., and Blankenship, D. D.: The 1500 m South Pole ice core: recovering a 40 ka environmental record, Ann. Glaciol., 55, 137–146,, 2014. 

Courville, Z. R., Albert, M. R., Fahnestock, M. A., Cathles, L. M., and Shuman, C. A.: Impacts of an accumulation hiatus on the physical properties of firn at a low-accumulation polar site, J. Geophys. Res.-Earth, 112, F02030,, 2007. 

Craig, H., Horibe, Y., and Sowers, T.: Gravitational Separation of Gases and Isotopes in Polar Ice Caps, Science, 242, 1675–1678,, 1988. 

Cuffey, K. M. and Paterson, W. S. B.: The Physics of Glaciers, Academic Press, 721 pp., 2010. 

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

Dahl-Jensen, D., Mosegaard, K., Gundestrup, N., Clow, G. D., Johnsen, S. J., Hansen, A. W., and Balling, N.: Past Temperatures Directly from the Greenland Ice Sheet, Science, 282, 268–271,, 1998. 

Dalrymple, P. C., Lettau, H. H., and Wollaston, S. H.: South Pole Micrometeorology Program: Data Analysis1, in: Studies in Antarctic Meteorology, American Geophysical Union (AGU), 13–57,, 1966. 

Denning, A. S., Fung, I. Y., and Randall, D.: Latitudinal gradient of atmospheric CO2 due to seasonal exchange with land biota, Nature, 376, 240–243,, 1995. 

De Rydt, J., Gudmundsson, G. H., Corr, H. F. J., and Christoffersen, P.: Surface undulations of Antarctic ice streams tightly controlled by bedrock topography, The Cryosphere, 7, 407–417,, 2013. 

Domine, F., Belke-Brea, M., Sarrazin, D., Arnaud, L., Barrere, M., and Poirier, M.: Soil moisture, wind speed and depth hoar formation in the Arctic snowpack, J. Glaciol., 64, 990–1002,, 2018. 

Dreyfus, G. B., Jouzel, J., Bender, M. L., Landais, A., Masson-Delmotte, V., and Leuenberger, M.: Firn processes and δ15N: potential for a gas-phase climate proxy, Quaternary Sci. Rev., 29, 28–42,, 2010. 

Epifanio, J. A., Brook, E. J., Buizert, C., Edwards, J. S., Sowers, T. A., Kahle, E. C., Severinghaus, J. P., Steig, E. J., Winski, D. A., Osterberg, E. C., Fudge, T. J., Aydin, M., Hood, E., Kalk, M., Kreutz, K. J., Ferris, D. G., and Kennedy, J. A.: The SP19 chronology for the South Pole Ice Core – Part 2: gas chronology, Δage, and smoothing of atmospheric records, Clim. Past, 16, 2431–2444,, 2020. 

Fahnestock, M. A., Shuman, C. A., Albert, M., and Scambos, T.: Satellite, Observational, Meteorological and Thermal Records from Two Sites in the Antarctic Megadunes – Stability of Atmospheric Forcing, Thermal Cracking, and the Seasonal Evolution of the Thermal Profile, AGU Fall Meeting Abstracts, C31C-03, 2004. 

Fudge, T. J., Biyani, S. C., Clemens-Sewall, D., and Hawley, R. L.: Constraining Geothermal Flux at Coastal Domes of the Ross Ice Sheet, Antarctica, 46, 13090–13098,, 2019. 

Fudge, T. J., Lilien, D. A., Koutnik, M., Conway, H., Stevens, C. M., Waddington, E. D., Steig, E. J., Schauer, A. J., and Holschuh, N.: Advection and non-climate impacts on the South Pole Ice Core, Clim. Past, 16, 819–832,, 2020. 

Giovinetto, M. B.: Glaciological Studies on the McMurdo-South Pole Traverse, 1960–1961, Research Foundation and the Institute of Polar Studies, The Ohio State University, 1963. 

Hamilton, G. S.: Topographic control of regional accumulation rate variability at South Pole and implications for ice-core interpretation, Ann. Glaciol., 39, 214–218,, 2004. 

Hawley, R. L., Courville, Z. R., Kehrl, L. M., Lutz, E. R., Osterberg, E. C., Overly, T. B., and Wong, G. J.: Recent accumulation variability in northwest Greenland from ground-penetrating radar and shallow cores along the Greenland Inland Traverse, J. Glaciol., 60, 375–382,, 2014. 

Herron, M. M. and Langway, C. C.: Firn Densification: An Empirical Model, J. Glaciol., 25, 373–385,, 1980. 

Huber, C., Leuenberger, M., Spahni, R., Flückiger, J., Schwander, J., Stocker, T. F., Johnsen, S., Landais, A., and Jouzel, J.: Isotope calibrated Greenland temperature record over Marine Isotope Stage 3 and its relation to CH4, Earth Planet. Sc. Lett., 243, 504–519,, 2006. 

Hudson, S. R. and Brandt, R. E.: A Look at the Surface-Based Temperature Inversion on the Antarctic Plateau, J. Climate, 18, 1673–1696,, 2005. 

Ikeda-Fukazawa, T., Hondoh, T., Fukumura, T., Fukazawa, H., and Mae, S.: Variation in N2/O2 ratio of occluded air in Dome Fuji antarctic ice, J. Geophys. Res.-Atmos., 106, 17799–17810,, 2001. 

Jordan, T. A., Martin, C., Ferraccioli, F., Matsuoka, K., Corr, H., Forsberg, R., Olesen, A., and Siegert, M.: Anomalously high geothermal flux near the South Pole, Sci. Rep., 8, 16785,, 2018. 

Jouzel, J., Barkov, N. I., Barnola, J. M., Bender, M., Chappellaz, J., Genthon, C., Kotlyakov, V. M., Lipenkov, V., Lorius, C., Petit, J. R., Raynaud, D., Raisbeck, G., Ritz, C., Sowers, T., Stievenard, M., Yiou, F., and Yiou, P.: Extending the Vostok ice-core record of palaeoclimate to the penultimate glacial period, Nature, 364, 407–412,, 1993. 

Jouzel, J., Masson-Delmotte, V., Cattani, O., Dreyfus, G., Falourd, S., Hoffmann, G., Minster, B., Nouet, J., Barnola, J. M., Chappellaz, J., Fischer, H., Gallet, J. C., Johnsen, S., Leuenberger, M., Loulergue, L., Luethi, D., Oerter, H., Parrenin, F., Raisbeck, G., Raynaud, D., Schilt, A., Schwander, J., Selmo, E., Souchez, R., Spahni, R., Stauffer, B., Steffensen, J. P., Stenni, B., Stocker, T. F., Tison, J. L., Werner, M., and Wolff, E. W.: Orbital and Millennial Antarctic Climate Variability over the Past 800,000 Years, Science, 317, 793–796,, 2007. 

Kahle, E. C., Steig, E. J., Jones, T. R., Fudge, T. J., Koutnik, M. R., Morris, V. A., Vaughn, B. H., Schauer, A. J., Stevens, C. M., Conway, H., Waddington, E. D., Buizert, C., Epifanio, J., and White, J. W. C.: Reconstruction of Temperature, Accumulation Rate, and Layer Thinning From an Ice Core at South Pole, Using a Statistical Inverse Method, J. Geophys. Res.-Atmos., 126, e2020JD033300,, 2021. 

Kawamura, K., Severinghaus, J. P., Ishidoya, S., Sugawara, S., Hashida, G., Motoyama, H., Fujii, Y., Aoki, S., and Nakazawa, T.: Convective mixing of air in firn at four polar sites, Earth Planet. Sc. Lett., 244, 672–682,, 2006. 

Kawamura, K., Severinghaus, J. P., Albert, M. R., Courville, Z. R., Fahnestock, M. A., Scambos, T., Shields, E., and Shuman, C. A.: Kinetic fractionation of gases by deep air convection in polar firn, Atmos. Chem. Phys., 13, 11141–11155,, 2013. 

Kindler, P., Guillevic, M., Baumgartner, M., Schwander, J., Landais, A., and Leuenberger, M.: Temperature reconstruction from 10 to 120 kyr b2k from the NGRIP ice core, Clim. Past, 10, 887–902,, 2014. 

Kobashi, T., Severinghaus, J. P., Brook, E. J., Barnola, J.-M., and Grachev, A. M.: Precise timing and characterization of abrupt climate change 8200 years ago from air trapped in polar ice, Quaternary Sci. Rev., 26, 1212–1222,, 2007. 

Kobashi, T., Severinghaus, J. P., and Kawamura, K.: Argon and nitrogen isotopes of trapped air in the GISP2 ice core during the Holocene epoch (0–11,500 B.P.): Methodology and implications for gas loss processes, Geochim. Cosmochim. Ac., 72, 4675–4686,, 2008. 

Kobashi, T., Kawamura, K., Severinghaus, J. P., Barnola, J.-M., Nakaegawa, T., Vinther, B. M., Johnsen, S. J., and Box, J. E.: High variability of Greenland surface temperature over the past 4000 years estimated from trapped air in an ice core, Geophys. Res. Lett., 38, L21501,, 2011. 

Landais, A., Barnola, J.-M., Masson-Delmotte, V., Jouzel, J., Chappellaz, J., Caillon, N., Huber, C., Leuenberger, M., and Johnsen, S. J.: A continuous record of temperature evolution over a sequence of Dansgaard-Oeschger events during Marine Isotopic Stage 4 (76 to 62 kyr BP), Geophys. Res. Lett., 31, L22211,, 2004. 

Landais, A., Barnola, J. M., Kawamura, K., Caillon, N., Delmotte, M., Van Ommen, T., Dreyfus, G., Jouzel, J., Masson-Delmotte, V., Minster, B., Freitag, J., Leuenberger, M., Schwander, J., Huber, C., Etheridge, D., and Morgan, V.: Firn-air δ15N in modern polar sites and glacial–interglacial ice: a model-data mismatch during glacial periods in Antarctica?, Quaternary Sci. Rev., 25, 49–62,, 2006. 

Leonard, K., Bell, R. E., Studinger, M., and Tremblay, B.: Anomalous accumulation rates in the Vostok ice-core resulting from ice flow over Lake Vostok, Geophys. Res. Lett., 31, L24401–L24401,, 2004. 

Leuenberger, M. and Lang, C.: Thermal diffusion: An important aspect in studies of static air columns such as firn air, sand dunes and soil air, nternational Atomic Energy Agency, Vienna, Austria, 2002. 

Lilien, D. A., Fudge, T. J., Koutnik, M. R., Conway, H., Osterberg, E. C., Ferris, D. G., Waddington, E. D., and Stevens, C. M.: Holocene Ice-Flow Speedup in the Vicinity of the South Pole, Geophys. Res. Lett., 45, 6557–6565,, 2018. 

Mariotti, A.: Atmospheric nitrogen is a reliable standard for natural 15 N abundance measurements, Nature, 303, 685–687,, 1983. 

Martinerie, P., Lipenkov, V. Y., Chappellaz, J., Barkov, N. I., and Lorius, C.: Air content paleo record in the Vostok ice core (Antarctica): A mixed record of climatic and glaciological parameters, J. Geophys. Res.-Atmos., 99, 10565–10576,, 1994. 

Miège, C., Forster, R. R., Box, J. E., Burgess, E. W., McConnell, J. R., Pasteris, D. R., and Spikes, V. B.: Southeast Greenland high accumulation rates derived from firn cores and ground-penetrating radar, 54, 322–332,, 2013. 

Morgan, J. and Severinghaus, J.: South Pole Ice Core Isotopes of N2 and Ar, USAP-DC [data set],, 2022. 

Orsi, A. J.: Temperature reconstruction at the West Antarctic Ice Sheet Divide, for the last millennium, from the combination of borehole temperature and inert gas isotope measurements, PhD Thesis, UC San Diego, 2013. 

Orsi, A. J., Cornuelle, B. D., and Severinghaus, J. P.: Magnitude and temporal evolution of Dansgaard–Oeschger event 8 abrupt temperature change inferred from nitrogen and argon isotopes in GISP2 ice using a new least-squares inversion, Earth Planet. Sc. Lett., 395, 81–90,, 2014. 

Oyabu, I., Kawamura, K., Uchida, T., Fujita, S., Kitamura, K., Hirabayashi, M., Aoki, S., Morimoto, S., Nakazawa, T., Severinghaus, J. P., and Morgan, J. D.: Fractionation of O2/N2 and Ar/N2 in the Antarctic ice sheet during bubble formation and bubble–clathrate hydrate transition from precise gas measurements of the Dome Fuji ice core, The Cryosphere, 15, 5529–5555,, 2021. 

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

Petrenko, V. V., Martinerie, P., Novelli, P., Etheridge, D. M., Levin, I., Wang, Z., Blunier, T., Chappellaz, J., Kaiser, J., Lang, P., Steele, L. P., Hammer, S., Mak, J., Langenfelds, R. L., Schwander, J., Severinghaus, J. P., Witrant, E., Petron, G., Battle, M. O., Forster, G., Sturges, W. T., Lamarque, J.-F., Steffen, K., and White, J. W. C.: A 60 yr record of atmospheric carbon monoxide reconstructed from Greenland firn air, Atmos. Chem. Phys., 13, 7567–7585,, 2013. 

Pietroni, I., Argentini, S., and Petenko, I.: One Year of Surface-Based Temperature Inversions at Dome C, Antarctica, Bound.-Lay. Meteorol., 150, 131–151,, 2014. 

Pollard, D., Chang, W., Haran, M., Applegate, P., and DeConto, R.: Large ensemble modeling of the last deglacial retreat of the West Antarctic Ice Sheet: comparison of simple and advanced statistical techniques, Geosci. Model Dev., 9, 1697–1723,, 2016. 

Price, P. B., Nagornov, O. V., Bay, R., Chirkin, D., He, Y., Miocinovic, P., Richards, A., Woschnagg, K., Koci, B., and Zagorodnov, V.: Temperature profile for glacial ice at the South Pole: Implications for life in a nearby subglacial lake, P. Natl. Acad. Sci. USA, 99, 7844–7847,, 2002. 

Schwander, J., Barnola, J.-M., Andrié, C., Leuenberger, M., Ludin, A., Raynaud, D., and Stauffer, B.: The age of the air in the firn and the ice at Summit, Greenland, J. Geophys. Res.-Atmos., 98, 2831–2838,, 1993. 

Severinghaus, J. P. and Battle, M. O.: Fractionation of gases in polar ice during bubble close-off: New constraints from firn air Ne, Kr and Xe observations, Earth Planet. Sc. Lett., 244, 474–500,, 2006. 

Severinghaus, J. P. and Brook, E. J.: Abrupt Climate Change at the End of the Last Glacial Period Inferred from Trapped Air in Polar Ice, 286, 930–934, 930–934,, 1999. 

Severinghaus, J. P., Sowers, T., Brook, E. J., Alley, R. B., and Bender, M. L.: Timing of abrupt climate change at the end of the Younger Dryas interval from thermally fractionated gases in polar ice, Nature, 391, 141–146,, 1998. 

Severinghaus, J. P., Grachev, A., and Battle, M.: Thermal fractionation of air in polar firn by seasonal temperature gradients, Geochem. Geophy. Geosy., 2, 1048,, 2001. 

Severinghaus, J. P., Grachev, A., Luz, B., and Caillon, N.: A method for precise measurement of argon 40/36 and krypton/argon ratios in trapped air in polar ice with applications to past firn thickness and abrupt climate change in Greenland and at Siple Dome, Antarctica, Geochim. Cosmochim. Ac., 67, 325–343,, 2003. 

Severinghaus, J. P., Albert, M. R., Courville, Z. R., Fahnestock, M. A., Kawamura, K., Montzka, S. A., Mühle, J., Scambos, T. A., Shields, E., Shuman, C. A., Suwa, M., Tans, P., and Weiss, R. F.: Deep air convection in the firn at a zero-accumulation site, central Antarctica, Earth Planet. Sc. Lett., 293, 359–367,, 2010. 

Shackleton, S., Bereiter, B., Baggenstos, D., Bauska, T. K., Brook, E. J., Marcott, S. A., and Severinghaus, J. P.: Is the Noble Gas-Based Rate of Ocean Warming During the Younger Dryas Overestimated?, Geophys. Res. Lett., 46, 5928–5936,, 2019. 

Souney, J. M., Twickler, M. S., Aydin, M., Steig, E. J., Fudge, T. J., Street, L. V., Nicewonger, M. R., Kahle, E. C., Johnson, J. A., Kuhl, T. W., Casey, K. A., Fegyveresi, J. M., Nunn, R. M., and Hargreaves, G. M.: Core handling, transportation and processing for the South Pole ice core (SPICEcore) project, Ann. Glaciol., 62, 118–130,, 2021. 

Sowers, T., Bender, M., Raynaud, D., and Korotkevich, Y. S.: δ15N of N2 in air trapped in polar ice: A tracer of gas transport in the firn and a possible constraint on ice age-gas age differences, J. Geophys. Res.-Atmos., 97, 15683–15697,, 1992. 

Sowers, T. A., Bender, M. L., and Raynaud, D.: Elemental and isotopic composition of occluded O2 and N2 in polar ice, J. Geophys. Res.-Atmos., 94, 5137–5150,, 1989. 

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

Sturm, M. and Johnson, J. B.: Natural convection in the subarctic snow cover, Quaternary Sci. Rev., 96, 11657–11671,, 1991. 

Suwa, M. and Bender, M. L.: Chronology of the Vostok ice core constrained by O2/N2 ratios of occluded air, and its implication for the Vostok climate records, Quaternary Sci. Rev., 27, 1093–1106,, 2008a. 

Suwa, M. and Bender, M. L.: O2/N2 ratios of occluded air in the GISP2 ice core, J. Geophys. Res.-Atmos., 113, D11119,, 2008b. 

Town, M. S., Waddington, E. D., Walden, V. P., and Warren, S. G.: Temperatures, heating rates and vapour pressures in near-surface snow at the South Pole, J. Glaciol., 54, 487–498,, 2008.  

Trudinger, C. M., Etheridge, D. M., Buizert, C., Hmiel, B., Krummel, P. B., Langenfelds, R. L., Manning, M., Mitrevski, B., Neff, P. D., Petrenko, V. V., Severinghaus, J. P., Smith, A. M., and Vollmer, M. K.: Rectifier Effect due to Seasonality in Convective Mixing in Firn, 2020, AGU Fall Meeting Abstracts, C033-02, 2020. 

Uemura, R., Motoyama, H., Masson-Delmotte, V., Jouzel, J., Kawamura, K., Goto-Azuma, K., Fujita, S., Kuramoto, T., Hirabayashi, M., Miyake, T., Ohno, H., Fujita, K., Abe-Ouchi, A., Iizuka, Y., Horikawa, S., Igarashi, M., Suzuki, K., Suzuki, T., and Fujii, Y.: Asynchrony between Antarctic temperature and CO2 associated with obliquity over the past 720,000 years, Nat. Commun., 9, 961,, 2018. 

van den Broeke, M. R. and van Lipzig, N. P. M.: Factors Controlling the Near-Surface Wind Field in Antarctica, Mon. Weather Rev., 131, 733–743,<0733:FCTNSW>2.0.CO;2, 2003. 

Van Liefferinge, B. and Pattyn, F.: Using ice-flow models to evaluate potential sites of million year-old ice in Antarctica, Clim. Past, 9, 2335–2345,, 2013. 

Verhulst, K. R.: Atmospheric Histories of Ethane and Carbon Monoxide from Polar Firn Air and Ice Cores, PhD Thesis, UC Irvine, 2014. 

Vihma, T., Tuovinen, E., and Savijärvi, H.: Interaction of katabatic winds and near-surface temperatures in the Antarctic, J. Geophys. Res.-Atmos., 116, D21119,, 2011. 

Waddington, E. D., Neumann, T. A., Koutnik, M. R., Marshall, H.-P., and Morse, D. L.: Inference of accumulation-rate patterns from deep layers in glaciers and ice sheets, J. Glaciol., 53, 694–712,, 2007. 

Weiler, K., Schwander, J., Leuenberger, M., Blunier, T., Mulvaney, R., Anderson, P. S., Salmon, R., and Sturges, W. T.: Seasonal Variations of Isotope Ratios and CO2 Concentrations in Firn Air, Science, 68, 247–272, 2009. 

Winski, D. A., Fudge, T. J., Ferris, D. G., Osterberg, E. C., Fegyveresi, J. M., Cole-Dai, J., Thundercloud, Z., Cox, T. S., Kreutz, K. J., Ortman, N., Buizert, C., Epifanio, J., Brook, E. J., Beaudette, R., Severinghaus, J., Sowers, T., Steig, E. J., Kahle, E. C., Jones, T. R., Morris, V., Aydin, M., Nicewonger, M. R., Casey, K. A., Alley, R. B., Waddington, E. D., Iverson, N. A., Dunbar, N. W., Bay, R. C., Souney, J. M., Sigl, M., and McConnell, J. R.: The SP19 chronology for the South Pole Ice Core – Part 1: volcanic matching and annual layer counting, Clim. Past, 15, 1793–1808,, 2019. 

Short summary
The composition of air bubbles in Antarctic ice cores records information about past changes in properties of the snowpack. We find that, near the South Pole, thinner snowpack in the past is often due to steeper surface topography, in which faster winds erode the snow and deposit it in flatter areas. The slope and wind seem to also cause a seasonal bias in the composition of air bubbles in the ice core. These findings will improve interpretation of other ice cores from places with steep slopes.